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ABSTRACT 

Ihis thesis investigates continuous-spatial dynamic models of 
the composition behavior of binary plate distillation columns. These 
partial differential equation models sire developed from basic distil- 
lation column principles and the discrete-plate models by treating 
the plate number as a continuous-spatial variable. Linearized 
continuous-spatial models are investigated in detail. 

The central purpose of this thesis is the development, pres- 
entation, suggested analytical solution technique, and example column 
evaluation of the Linear Polynomial -Coefficient Model (LPCM) which 
is a linearized continuous-spatial model in which the coefficients 
in the partial differential equation are n-th degree polynomials in 
the spatial variable. A general analytical solution technique is 
proposed in which the spatial differential eigenvalue problem 
resulting from separation of variables is transformed to a Liouville 
Normal-Form equation which is then converted to a homogeneous 
Fredholm II integral equation. Several simple examples of the model 
are solved in complete detail, and the proposed solution technique 
is applied to a model with first-degree polynomial coefficients. 

An analytical and a computational analysis giving the details of 
each step of the solution technique is presented. It is suggested 
that greatly reduced computation times compared to discrete models 
will result from application of the proposed solution technique, 

A bibliography of 352 references, 202 of which pertain directly 
to distillation column dynamics and control, is presented and related 
to the areas of the thesis, 
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SECTION 1 



INTRODUCTION TO DISTILLATION (i) 

11 BASIC PRINCIPLES OF BINARY DISTILLATION 

12 THE BINARY PLATE DISTILLATION COLUMN 

13 PHILOSOPHY AND DISTILLATION 

14 A BRIEF REVIEW OF THE LITERATURE OF DISTILLATION 

THE FUNDAMENTAL AND UNIVERSAL ABSOLUTE: 

"EXISTENCE, REALITY, THE EXTERNAL WORLD, IS WHAT IT IS, INDEPEN- 
DENT OF MAN’S CONSCIOUSNESS, INDEPENDENT OF ANYONE'S KNOWLEDGE, 

JUDGMENT, BELIEFS, HOPES, WISHES, OR FEARS - THAT FACTS ARE FACTS, 

THAT A IS A, -THAT THINGS ARE WHAT THEY ARE." 

NATHANIEL BRANDEN (B-28) 

THIS SECTION PRESENTS THE BASIC PHYSICAL AND MATHEMATICAL IDEAS 
OF BINARY DISTILLATION ORIENTED TOWARD FORMING A MATHEMATICAL MODEL 
OF A BINARY PLATE DISTILLATION COLUMN. THE GENERAL PROBLEM OF OBSERV- 
ING AND MODELING PHYSICAL SYSTEMS IS DISCUSSED AND THE HISTORY AND 
GENERAL LITERATURE OF DISTILLATION ARE REVIEWED BRIEFLY. SEVERAL 
OPINIONS ON THE PHILOSOPHY OF DESERVING REALITY AS RELATED TO MODELING 
A DISTILLATION COLUMN AND ON THE DESIRE TO ACHIEVE PROFIT AS RELATED 
TO THE COST OF SEPARATION OF COMPONENTS ARE PRESENTED. 

THE CENTRAL PURPOSE OF THIS THESIS IS THE DEVELOPMENT, PRESENTATION, 
SUGGESTED SOLUTION TECHNIQUE, AND EVALUATION OF THE LINEAR POLYNOMIAL- 
COEFFICIENT MODEL (LPCM) OF THE DYNAMIC BEHAVIOR OF A BINARY PLATE 
DISTILLATION COLUMN, THIS SECTION PRESENTS THE BASIC PRINCIPLES LEADING 
UP TO THE COMPLETE DEVELOPMENT OF THE LPCM IN SECTION 2(M) AND THE 

PRESENTATION OF AN INTEGRAL EQUATION SOLUTION TECHNIQUE IN SECTION3(L). 
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CHAPTER II 



BASIC PRINCIPLES OF BINARY DISTILLATION 

The purpose of Chapters II and 12 is to present the basic 
principles behind the separation of binary mixtures and to develop 
the techniques for describing mathematically and graphically the 
steady state characteristics of a distillation column. These topics 
are presented for several reasons: 

1. Some readers may not be acquainted with the basic theory, 
which is necessary for the developments which follow. 

2. This is a convenient method for developing a consistent 
notation for use throughout this thesis. 

3. Presentations of those basic principles specifically 
used in this thesis are not often found in the existing literature 
in forms easily understood by readers without some background in 
the subject. 

There are many available techniques for describing binary 
distillation in quantitative or qualitative terms. Most of these 
techniques are oriented specifically toward one of these two outlooks. 
The technique to be described in Chapters II and 12 is the use of the 
McCabe-Thiele diagram which has the unique advantage of a quantitative, 
qualitative, and, in a sense, visual insight into binary distillation. 

The theory developed in this chapter represents a combination 
of the developments available in the literature. The textbooks used 
in this presentation of the theory in this chapter are listed below 
in the order of decreasing pertinence to this chapter. 

Gould (G-3) 
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Holland 



Bennett and Myers 



Van Winkle 



(B-l) 

(V-l) 

(H-2),(H-7) 



Some of the other textbooks which present the basic principles of 
distillation are: (T-l), (S-l), (H-l), (C-l), (R-17), (C-5), (H-6), 

(H-5), (R-21), and (M-2). 

II. 1 DEFINITION OF BINARY DISTILLATION AND RELATIVE VOLATILITY 

Binary distillation is defined as a process which separates a 
mixture of two components by utilizing mass transfer between the 
liquid and vapor phases of the components. The essence of this sep- 
aration lies in the fact that when the vapor and liquid phases of a 
binary mixture are in equilibrium, the vapor is richer in the lighter 
component than is the liquid. The process by which the vapor phase 
becomes richer in the component which boils at the lower temperature 
(lighter component) is called mass transfer. Although the term 
distillation is occasionally used to describe the removal of volatile 
materials from solids, the term as used in this thesis will apply only 
to the separation of volatile components found in liquid solutions. 

The equilibrium mentioned above is a phase equilibrium in which 
the properties of the two phases depend upon the physical characteristics 
of the two components which are present. The primary physical char- 
acteristic of interest in this thesis is the relative volatility a* 

If the concentration of the lighter component in the liquid phase is 
defined as u and the concentration of the lighter component in the 
vapor phase is defined as f(u), then the relative volatility is given 
by equation II. 1. This expression is only valid when a does not vary 
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a " izH . f Cu) 
u i-f(u) 



II. 1 



with composition, which is a valid approximation for a large number 
of liquids. Equation II. 1, then, defines equilibrium between the two 
phases in terms of the constant relative volatility a. 

11.2 THE EQUILIBRIUM CURVE 

The assumption that a is constant is essential to the develop- 
ments of this thesis because it allows a specific equilibrium function 
to represent the characteristics of the phase mixture. Thus, an 
equilibrium phase diagram can be drawn as in Figure 11,1 and a specific 
equilibrium function f(u) can be expressed in Equation 11.2. 

f(u) - au 11.2 

~Ta-Du 




Figure 11,1 - THE EQUILIBRIUM CURVE 
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In Figure II. 1 the curved line (equilibrium curve) represents 
a graphical plot of equation 11.2. If the relative volatility is not 
assumed to be constant, then this curve may have a significantly 
different shape and be represented by a different functional relation- 
ship. In trying to understand the meaning of the equilibrium curve 
it is helpful to consider a point u Q in Figure 11,1, 

Given a binary mixture which has liquid and vapor phases at 
equilibrium, then the concentration of the lighter component in the 
liquid is given by u Q , and the concentration of the lighter component 
in the vapor is given by f(u Q ), Thus, going from the liquid to the 
vapor in the two phase mixture corresponds to going from the u line 
to the f(u) line on the equilibrium curve. The fact that the relative 
volatility is such that f(u) > u in this case implies that a partial 
separation of the mixture can be accomplished by separating the vapor 
from the liquid. 

11.3 USING THE EQUILIBRIUM CURVE TO DESCRIBE SEPARATION 

If the vapor is separated from the liquid in a two phase binary 
mixture, then condensing the vapor produces two separate liquids of 
different compositions, the one having been condensed being richer in 
the lighter component. This method of separation is the key to binary 
distillation. 

This separation method can be visualized by using the equilibrium 
curve in Figure 11.2, The original liquid at composition u q is boiled 
and part of it becomes vapor at composition f(u Q ). If the vapor at 
f(u Q ) is then separated from the liquid at u q and then condensed, the 
liquid condensate is of composition u^ = f(u Q ). The new liquid at 
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composition u-^ can now be boiled to produce vapor of composition 
f(u^), where f(u^) > u^. Thus, by successively boiling and condensing 
the liquids and vapors, a desired degree of separation in terms of 
the upper and lower (on the graph of Figure 11,2) liquids can be achieved. 
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11.4 PHYSICAL SEPARATION OF VAPOR AND LI CHID 



The next problem to be considered is that of the physical equip- 
ment necessary to separate the vapor of f(u Q ) from the liquid at u , 
while at the same time allowing the liquid-vapor contact between the 
vapor at f(u Q ) and the liquid at u^. There are two commonly used 
types of equipment for accomplishing this: packed columns and plate 

columns, both of which have innumerable variations. 

Packed columns use a form of continuous contacting of liquid 
and vapor phases. This is accomplished by use of many small devices 
in the form of rings, saddles, spheres, etc, , which are randomly 
"packed" into the column. This arrangement is designed to provide a 
very large surface area for a given tower volume because the contacting 
and mass transfer occur at the surfaces of the particles. The 
mathematical models (Section M) of packed columns are usually partial- 
differential equations. Packed columns will not be investigated 
specifically in this thesis, but it Is expected that the developments 
of this thesis could be applied to packed columns, 

Plate columns accomplish vapor-liquid separation and vapor-liquid 
contacting by allowing the liquid at composition u^ to flow over the 
top of a plate, which is designed to let the vapor at f(u Q ) pass 
through into the liquid flowing over it. There are a multitude of 
different plate designs which accomplish this efficiently, such as 
bubble-cap, perforated, and sieve plates. Figure 11,3 shows the basic 
physical characteristics of a bubble-cap plate, Ihe concentrations 
shown in Figure 11.3 correspond to those on the equilibrium curve of 
Figure 11.2, Using these two figures one can see both physically and 
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f 




Figure 11.3 - OPERATION OF A BUBBLE-CAP PLATE 



conceptually how separation is achieved. 

Physically, the liquid at u-^ enters the plate of Figure 11.3 from 
the left-side downcomer of the plate above. The u^ liquid then flows 
across the plate, in this case left to right, and is mixed (contacted) 
with the vapor f(u Q ) coming from beneath the plate. These two phases 
reach an equilibrium such that the liquid leaving the plate is at u Q 
and the vapor leaving the liquid is at f(u^). Since f(u^) > u^, a 
partial separation of the components results. Conceptually, following 
these same arguments on the lines in Figure 11,2 reveals the properties 

of the mechanism of separation in terms of the mathematical expressions 
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describing equilibrium. 

This chapter has presented the basic principles of binary 
separation using distillation. This separation has been described 
conceptually in terms of the relative volatility and the resulting 
equilibrium curve and physically in terms of the operation of a 
bubble-cap plate. The next chapter puts several of these plates 
together and presents the physical and conceptual descriptions of 
a complete distillation column. 
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CHAPTER 12 



THE BINARY PLATE DISTILLATION COLUMN 

The particular binary separation device chosen for study in this 
thesis is the plate distillation column. A plate distillation column 
is a vertical cascaded arrangement of individual plates (Figure 11.3) 
supported in a tower or column. The distillation column alone cannot 
perform separation but requires auxiliary equipment for its operation. 
This chapter presents a very brief description of the physical operation 
of a distillation column and its auxiliary equipment, a graphical 
description of its operation in terms of the McCabe-Thiele diagram, 
and a steady-state mathematical model of its operation. 

12.1 PHYSICAL OPERATION OF A BINARY PLATE DISTILLATION COLUMN 

The physical separation of a binary mixture by an individual plate 
was described in Chapter II. In general, one plate is usually in- 
adequate to achieve the desired degree of separation of the mixture, 
thus many plates are combined in cascade to achieve greater output 
purity. An arrangement of eleven such plates is shown in the column 
of Figure 12,1. On any given plate in this column the liquid and vapor 
mix and reach equilibrium, with the vapor rising in the column plate- 
by-plate through the bubble caps and with the liquid flowing back and 
forth down the column. 

Two of the auxiliary equipments necessary to the operation of the 
column are the condenser and the reboiler. The vapor flowing out of 
the top of the column enters the condenser where it is liquified and the 
resulting liquid is partially fed back to the top tray and partially 
removed as tops product. Similarly, the liquid flowing out the bottom 
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of the column is divided into that portion removed sis bottoms product 
sind that portion which is vaporized in the reboiler and fed back to the 
bottom plate. The other auxiliary equipments, such as pumps to cir- 
culate the fluid, valves for control, structural supports, and many 
others, are not shown in Figure 12.1 but are, nevertheless, essential 
to the operation of the system. 

Having described the operation of the internal cycling of the 
fluids in the column, the overall operation of the system can now be 
discussed. Ihe input mixture to be separated is fed to that tray 
designated as the feed tray and enters the fluid cycle. The separated 
outputs are then taJcen off as the top and bottom products. Energy for 
operation of the system is supplied by the reboiler, steam heated in 
the case of the column of Figure 12.1, and energy is removed from the 
system by the condenser, ' water cooled in Figure 12.1. 

There are innumerable types, arrangements, designs, and sizes of 
distillation columns, just as there are many different aspects of any 
given column which can be studied, such as chemical, thermal, structural, 
fluidic, economic, and environmental aspects. Ihe central purpose of 
any column is to separate a binary mixture and the main criterion of 
"goodness" of the operation of the column is how well it performs this 
separation. Distillation columns are usually designed to separate the 
mixture to desired purity in such a way as to maximize the economic 
profit derived from the sale of the separated products. 

Output variations in the product purity greatly affect the profits 
derived from such sales. If the output is overly pure, then the revenue 
resulting from the sale of material thought to be of lower purity will 
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not be as high as it could have been. If the output is lower than the 
desired composition, then it cannot be sold at the price set for 
material of higher quality and profits will be reduced, Thus, composi- 
tion variations in distillation column inputs and outputs are rather 
important, and mathematical models for the interaction between feed 
composition changes and output composition changes will be developed 
in this thesis for the purpose of studying such variations. 

12.2 THE MCCABE- THIELE DIAGRAM 

One of the best graphical techniques for describing the steady- 
state operation of a binary plate column is the McCabe-Thiele diagram 
(M-10). A technique for describing the steady-state operation of the 
column is important to considerations involving variations in composi- 
tions, especially when composition variations are to be interpreted 
as upsets or deviations from an Initial steady state to a final steady 
state. The basis of the McCabe-Thiele diagram is the equilibrium curve 
presented in Figure 11.2 for the operation of one plate. For the eleven 
(ll) plate column of Figure 12.1 the McCabe-Thiele diagram for a = 3»0 
(M-15) is shown in Figure 12.2. 

The graphical expression of the individual plate compositions 
presented by the McCabe-Thiele diagram gives both a qualitative and a 
quantitative view of the steady-state operation of the column. Qualita- 
tively, it shows how the individual plates act to increase the purities 
of the output streams by operating in cascade with each step up the 
diagram representing a plate corresponding to a step up the actual 
column. Quantitatively, the individual plate steady-state compositions 
and the operating line slopes can be read directly from the diagram. 
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Figure 12.2 

MCCABE-THIELE DIAGRAM FOR 
AN ELEVEN PLATE COLUMN (M-15) 
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The actual operating characteristics of the column are evident 
from the McCabe-Thiele Diagram. The amounts or percentages of liquid 
and vapor reflux flows determine the slopes of the two operating lines, 
the upper and lower lines in Figure 12.2. Each point of intersection 
on the operating lines represents the liquid composition of a plate, 
and thus the number of liquid points equals the number of plates in 
the column. Each point of intersection on the equilibrium curve repre- 
sents the vapor composition between two of the column plates. The 
intersection of the two operating lines is determined by the q-line 
which depends upon the properties and condition of the feed at u^. All 
of these properties are represented by quantities and equations in the 
steady-state model of the distillation column, 

12.3 A STEADY-STATE MODEL OF A BINARY PLATE DISTILLATION COLUMN 

This subsection begins the developments leading to mathematical 
equations describing the operation of a binary plate distillation column. 
A mathematical model of the steady-state operation is to be presented 
and explained. Since most of the detailed developments of steady-state 
models are well covered in the literature (See B2.1 - Steady-State 
Analysis and McCabe-Thiele Diagrams), the discrete-plate steady-state 
model and its symbols will merely be presented, not derived in detail, 
in this section. The format of this presentation will be to present 
the description of the symbols to be used, present the model, and give 
a brief explanation of the model and the assumptions behind it. 

The presentation of the discrete-plate steady-state model begins 
by listing and describing the mathematical symbols to be used in the 
model. This list of symbols is found in Table 12.1. The symbols are 
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12.1 



E n “ f ' K) ~ u n 

f( u n) - 

Where f ' (u n ) is a non-equilibrium state 

u - Concentration of the lighter component in the liquid (lbm 
mole lighter comp. /lbm liquid) 

f(u) - Concentration of the lighter component in the vapor (lbm 
mole lighter comp. /lbm vapor) 

a - Relative volatility (dimensionless); assumed constant 

F - Feed rate (lbm liquid/hour) 

D - Distillate rate (lbm liquid/hour) 

W - Withdrawal rate of Bottoms Product (lbm liquid/hour) 

,L - Liquid rate in the upper section (lbm liquid/hour); assumed 
constant 

- Liquid rate in the lower section (lbm liquid/hour) ; assumed 
constant 

V - Vapor rate in the column (lbm vapor/hour) ; assumed constant 
B u - Upper reflux ratio (lbm liquid/lbm vapor), L u /V 

- Lower reflux ratio (lbm liquid/lbm vapor), L^/V 
k - Feed plate index (number) ; integer 

n - Internal plate index (number); integer 1 < n < N 

q - Portion of the feed which adds to the lower liquid rate; 

L, « L + qF 
1 u u 

Table 12.1 - LIST OF SYMBOLS USED IN DISTILLATION 
COLUMN STEADY-STATE MODELS 
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fairly standard throughout the literature with one important exception. 
Usually, in the literature, the composition of the lighter component 
in the liquid is represented by x and that in the vapor by y. In this 
thesis x will be used in Sections M and L to represent a continuous- 
spatial variable or continuous-plate-number variable; and y is not used, 
but the term f(u) is used instead. The term f(u) is standard in most 
of the literature. The symbols in Table 12.1 also apply to the corres- 
ponding symbols placed on Figure 12.1. 

Figure 12.3 presents the discrete-plate steady-state model of a 
binary plate distillation column. This model is merely a combination 
of the conservation-of-mass equations (inflow = Outflow) for each plate 
and a listing of the end, feed, and equilibrium description. These 
equations say mathematically what the McCabe-Thiele diagram says graph- 
ically. Solving these equations mathematically for the plate composi- 
tions is also equivalent to "stepping off" the compositions on the 
McCabe-Thiele diagram. Tne amount of mathematical manipulation involved 
in solving for the plate compositions when there are many trays is 
rather large when compared with the ease of stepping off these numbers 
on the McCabe-Thiele diagram. 

Several very significant assumptions are implicit in the model 
statement of Figure 12.3. The upper and lower liquid rates, the 
column vapor rate, and the relative volatility have all been assumed 
to remain constant. In addition, the Murphree (M-8) efficiencies 
defined in equation 12.1 have all been assumed to be unity, or equiv- 
alently, the assumption is made that each plate operates at complete 
equilibrium, 
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Model Equations 



Rectification Section (Upper) 

f <Vl> ■ Vn + 

Stripping Section (Lower) 

f(u n-l> * Vn + ( 1 “ B l> u w 

Where B = L /V and B. = L /V 
u u' 1 i' 



k + 1 < n £ N 



1 < n ^ k 



Equilibrium 



f(0 



aU n 



n l+(a-l)u n 



End Conditions 



f( u Ij) “ u d 


Top 


n 


- N 




" "f 


Feed 


n 


- Ji. 


U 1 


= “v 


Bottom 


n 


= 1 



Feed Condition 

= L u + qF or B^ = B u + qF/V ; q = 1 in this case. 
Figure 12.3 - DISCRETE-PLATE STEADY-STATE MODEL 



To summarize briefly, Chapters II and 12 present the basic 
principles and mathematical descriptions for the operation of a binary 
plate distillation column. This presentation takes the form of a 
graphical description in terms of equilibrium curves and the McCabe- 
Thiele diagram, a brief physical description, and a set of mathematical 

equations forming a discrete-plate model for the steady-state operation. 
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The mathematical developments of this thesis in terms of distilla- 
tion column dynamics or transient models continues in Section M. In 
that section, the discrete-plate steady-state model in Figure 12.3 is 
seen to be a subcase of the discrete-plate dynamic model, and a continu- 
ous-spatial steady-state model (not presented) is seen to be a subcase 
of a continuous-spatial dynamic model. 
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CHAPTER 13 



PHILOSOPHY AND DISTILLATION 

This chapter is purely opinion and will be written using the 
first person. It seems to me that some comments can be made about 
two very general philosophical desires as applied to the process of 
distillation. Ihe first of these is the desire to observe and explain 
reality and the second is the desire to achieve profit. Philosophy 
is the study of the principles of reality, and distillation is a 
process for separating components; these two would seem to be almost 
completely unrelated, but, on the contrary, I think a definite and 
crucial relationship exists and propose to present it in this chapter. 

13.1 EXISTENCE EXISTS ! 

The fundamental axiom of reality is that existence exists, or as 
Aristotle stated, "A is A” (B-28). I exist in this reality as a being 
of volitional consciousness with only two choices open to me: to live 

or to die. I choose to live. I cannot live except by acquiring knowledge 
of one form or another. My only means of acquiring knowledge is through 
my senses. The only tool which I have for the thinking required to 
acquire knowledge from the inputs of my senses is the ability to reason. 
Reason requires the use of logic which is "the art of non-contradictory 
identification" (Atlas Shrugged - Ayn Rand). Thus, I wish to acquire 
knowledge about reality by observing it and applying logic to those 
observations. 

One of the many questions I might ask is, "Does there exist a 
defined limit to the amount of knowledge attainable from the observations 
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of my senses and any devices which I may design to aid them?” If 
I wish to observe the operation of a distillation column, can I say 
beforehand that I'll never be able to understand completely the exact 
mechanisms or motions undergone to achieve the desired separation of 
components? I think that there are several practical resaons why 
the answers to these questions for an individual must be YES, even 
though the philosophical axiom A = A implies that the answers must 
be NO. 

The first practical constraint that I run into is my limited 
physical capacity for knowledge in terms of my limited lifetime for 
acquiring it. But suppose for a moment that I am granted an infinite 
lifespan in which to observe and apply logic to distillation columns. 

I might then attempt to make more specific models of the column based 
upon more accurate observations of its operation. At some point, 
however, I run into another practical constraint: my act of observing 

the operation of the column begins to affect the operation which I'm 
trying to understand. This would be some form of "uncertainty principle" 
applied to distillation column measurements. 

Suppose I don't accept uncertainty; after all, I'm certain that 
A = A. This supposition means that I can eventually determine a model 
of the operation of a distillation column that tells me to N -» oo decimal 
places exactly what each molecule, atom, or even sub- atomic particle 
is doing within the column at any time, including the effects of my 
observations, Tnus, philosophically, the axiom A = A implies that 
unlimited reason is competent to define reality. 
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The certainty of limited life spans of individuals imposes a 
practical constraint on the amount of knowledge about distillation 
attainable by any one man. The assumption that human life will 
propagate forever through time implies that unlimited reason exists. 
Accepting this assumption one can say that man, now and future, has 
the unlimited reason necessary to define reality. What, then, deter- 
mines the portion of the limited facilities of an individual which will 
be applied to the study of distillation? This, then, gets into the 
area of the second philosophical desire, or necessity: to achieve profit 

in order to live, 

13.2 INDIVIDUAL PROFIT MAXIMIZATION 

The amount of effort devoted by any person to the study of distil- 
lation will be in proportion to the benefit (or profit) returned to 
that individual as a result of his efforts. For some this benefit 
may be derived by placing high value upon the self-satisfaction of 
explaining, even partially, the operation of a complex physical system, 
i.e. knowledge for the value of knowledge. Knowledge contains no life 
sustinence, and one who places a very large personal value in knowledge 
must have other means of attaining that profit convertible to life 
sustinence, either by selling that knowledge for food or by devoting 
a portion of his efforts, which could have been spent on studying 
distillation, to the acquisition of food by some other endeavor. 

It is my opinion that each individual attempts to maximize the 
value to him (profit) of the use of his intellect. The individual who 
receives a small (but adequate for life support) income and devotes 
the remainder of his life to watching television is maximizing his 
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value implicitly by placing high value on the endeavor of watching 
television. However, in order for such an individual to have received 
that income necessary for life support, he must have devoted energies 
to some endeavor in which the profits are transformable into food. 

Thus, the conclusion can be made that no matter how an individual 
seeks to maximize the returns from his endeavors some portion of his 
life must be devoted to endeavors which produce something which can be 
transformed into food. 

The question of whether or not there will definitely be people 
throughout time who will devote a portion of their life to studying 
distillation then becomes equivalent to asking whether or not distil- 
lation is an endeavor which produces something which can be transformed 
into food. If distillation is a profitable endeavor, then the relative 
question of how much of an individual's life will be spent studying it 
is answered by defining just how profitable distillation is to the man 
who studies it. 

13.3 PROFIT IN CERTAINTY 

A distillation column has an input, which is a mixture of two com- 
ponents, and two outputs each of which, one is certain, contains a 
higher concentration of a given component of the mixture. It is a 
fact that in the present world the total value of the two outputs is 
greater than the sum of the mixture value and the value given up to 
perform the separation. Thus, distillation IS a profitable endeavor. 
WHY? 

Distillation is profitable because certainty is of very high 
value. The operation of a distillation column produces two outputs 



- 35 - 



of greater certainty than the input. The fact that the composition 
of the outputs is certain is directly transformable into values 
equivalent to food and, as such, is a profitable endeavor for an 
individual to use for life support. The existence of profit from 
some endeavor implies that the benefits derived from it are greater 
than the costs incurred. 

What is the cost of certainty? Suppose that I knew a microscopic 
man capable of distinguishing between the two different molecules of 
a binary mixture and capable of deflecting one type of molecule in one 
direction and the other type in another direction. My microscopic 
friend could then be the Maxwell’s Demon of distillation if I let him 
stand in the inlet of the feed pipe to the distillation column and 
bat the lighter molecules upward and let the heavier ones fall, sepa- 
rating the mixture for me. I could then shut down the reboiler and 
condenser and sell the resulting pure outputs with no cost to me in 
terms of energy supplied to the distillation column. In fact, I could 
reason that the only effort involved in the entire process is the intel- 
lectual effort exerted by my small friend in recognizing which molecules 
are which. Thus, the cost of certainty is the information supplied by 
my friend, 

A supposed fallacy of this argument is often presented as follows. 
The second law of thermodynamics states that entropy is always increas- 
ing in a real, physical process; entropy is then equated to average 
uncertainty, and uncertainty is defined as loss of information. So, 
somewhere down the line, someone loses out; in this case it was, 
according to this argument, my small friend. Ihe Demon's information 
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was being transformed into my profit. The cost of certainty is 
increased uncertainty, if one accepts the concept of average uncer- 
tainty. 

Suppose, once again, that I state that my concept of existence 
denies uncertainty the right to exist for all time. The concept NOT 
is the negation of reality and, as such, cannot exist. No-one will 
ever observe the color Not-Blue because it cannot exist as such. 

Ihe use of the concept of average uncertainty is a denial of the 
statement A = A and can only be rejected as false if A = A is accepted 
as true. The average uncertainty argument implies that uncertainty is 
always increasing; A = A implies that certainty is increasing. The 
two are diametrically opposed, I accept A = A and deny that the 
universe is running down. 

What then is the cost of the separation achieved by my small 
friend if I reject the concept of average uncertainty? The cost is 
the certainty used by my friend in separating the mixture. The profit 
which I derive from the sale of the two separated components results 
from the fact that the two separated components have greater certainty 
and therefore greater value than the value of the certainty used by 
ray friend in separating the mixture. Thus, I can pay my friend for 
his services and we'll both be better off as a result of the process. 
Thus, all profit is the result of certainty. 

13.4 CERTAINTY AND KNOWLEDGE 

Ihe only route to certainty is through knowledge. The only route 
to knowledge is through observation and reason utilizing logic. If 
I desire certainty about the operation of a distillation column, then 
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I must turn one on (or read about someone who did), watch it, and apply 
reason to my observations to acquire knowledge about it. My limited 
physical capacity and lifespan dictates that my optimum acquisition of 
knowledge about distillation occurs during the time span up to the 
point where my marginal acquisition of certainty for expenditure of 
concentration begins to become negative. Beyond that point I would be 
better off studying some other profitable venture. Thus, constraints 
must be imposed on the acquisition of knowledge about distillation. 

The physical capacity constraints of an individual are realisti- 
cally applied to the study of distillation through judicious approxi- 
mation. I realize that given enough time and intellectual capacity 
I could explain distillation to any desired completeness. However, in 
order to acquire any knowledge about distillation to use in deriving 
profit from my study, without already understanding it completely, I 
must make simplifying approximations and build up a hierarchy of 
knowledge pertaining to the subject. This is exactly what is done, 
and the following sections in this thesis present a theory which over- 
flows with simplifications and approximations. The final test of any 
theory based on approximations must be that the revenue resulting 
from application of the theory is greater than the expenditure of 
effort in developing it. Such a test has yet to be applied to the 
theories presented in this thesis. 

13.5 WHAT'D HE SAY? 

This chapter has presented my opinions concerning the relation- 
ship between my concept of philosophy and the study of distillation. 
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The chapter begins with arguments which position me as a human being 
within the framework of reality as a being who must know to exist. 

The reason I must know is that I must acquire profit (food) to live. 
Distillation is shown to be profitable by arguments which show that 
distillation increases certainty (value) and by arguments which show 
that uncertainty cannot exist if one assumes that existence exists. 
The same arguments are used to invalidate the uncertainty principle 
and the concept of average uncertainty. Finally, the use of simplifi 
cations and approximations in the study of distillation is justified 
on the basis that a man has limited facilities for use in studying it 
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CHAPTER 14 



A BRIEF REVIEW OF THE LITERATURE OF DISTILLATION 

Hie volume of literature available pertaining to the study of 
distillation is completely overwhelming! It is easy to see why there 
can be no "Renaissance Man" in modern times. Hie bibliography pre- 
sented in Chapter B1 of this thesis contains 352 references of which 
202 pertain to distillation column dynamics and control. Hie 352 
references are minute in number compared to the total literature. 

The 202 references represent a significant percentage of the total 
material available pertaining to column dynamics and control. By 
far the greatest proportion of the total literature of distillation 
deals with column operation, design, and steady-state analysis. 

Distillation is such a general area and has so many inter- 
actions with other areas that it is easy to see why so much can be 
said about it. Hie areas listed in Chapter B2 represent, in themselves, 
an extensive collection of knowledge, and the list is fax from complete. 
In this thesis the term "literature" will refer to this infinitisimal 
subset of 352 references. 

Hie purpose of this chapter is to present descriptions, comments, 
and notes pertaining to some of the references in Chapter Bl. With 
the exception of some of the theses, each entry in the bibliography 
was very cursorily inspected by the author, and some notes and comments 
were made where deemed worthwhile. This chapter is a connected 
listing of those notes presented in the format of Table B2.1, No 
attempt has been made to list complete descriptions, and the author 
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has chosen to use abbreviated language (incomplete sentences, minimum 
articles, etc.) to increase the densite of information with, hopefully, 
minimum loss of comprehension. 

Many of these comments represent this author's opinions and should 

be evaluated as such by the reader. Those references in Chapter B2 

under each area about which no information is given were most likely i 

a. Dynamic response or control theses not read (not available) 

b. Poorly written or invalid (author's opinion) 

c. Not understood by the author after reading, but seemingly 
applicable 

d. Those presenting material covered better in other references 
(author's opinion) 

e. Partially explained by the title 

f. Not really significant to the area but containing pertinent 
material 

g. Discussed under another area in which case the area will be 
listed. 

The last name of the first author of each reference is listed for the 
mnemonic convenience of the reader. In the interest of space conserva- 
tion without loss of meaning all area titles and subtitles in Table 
B2.1 have been numbered and brought to the margin. 

14.1 GENERAL THEORY OF DISTILLATION 
1. Textbooks 

(B-l) Bennett & Meyers - covers heat, mass, momentum transfer, 

binary and multicomponent distillation, stagewise operations - 
primarily chemical engineering presentation; good presentation 
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of basic principles of distillation and steady-state analysis. 

(G-3) Gould - very general process control text; presents detailed 
mathematical methods, models, and analysis for the dynamic 
behavior and control of a large number of chemical processes, 
including packed and plate columns. 

(H-10) Henley & Staff in - very clear presentation of basic principles 
of distillation. 

(R-23) Reid & Sherwood - presents in great detail the equilibrium, 
diffusion, thermodynamic , and thermal properties of liquids 
and gases. 

(R-2l) Robinson & Gilliland - detailed presentation of basic principles 
and steady-state analysis. 

(V-2) Van Winkle - guide to fractional distillation design; considers 
hydraulics, multicomponent aspects, all sites of distillation 
design of plate and packed columns; discusses history of distil- 
lation; see also 3» 

(A-19) Aris - steady-state distillation analysis. 

(B-15) Bodman - presents Fortran IV program for optimum design of 

stagewise - ethylbenzene vacuum distillation reactor; economic 
optimization. 

(M-14) McCabe & Smith - recent text presenting basic principles, 

(0-2) Oliver - recent text presenting basic principles. 

(H-2) Holland - uses Ihiele-Geddes calculational procedure and Q 

method of convergence; mostly steady-state; "Instead of seeking 
exact analytical solutions for models that roughly approximate 
the actual system, researchers have put a vast amount of effort 
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into the development of iterative procedures in which 
progressively better initial values of the independent variables 



are selected for each successive trial." 

(H-9) - 8. 

(H-18) Hengstebeck - Considers many aspects of column design; pre- 
sents basic principles and steady-state analysis; practical 
text. 

(L-26) - 9. 

(C-l) Campbell - steady-state analysis; dynamic analysis using 
signal flow graphs, frequency analysis; process control. 

(P-3) Peters - plant design using detailed plate and bubble cap 

design equations; economics of distillation; optimum reflux 
design. 

(P-10) Pratt - recent text presenting basic principles, steady-state 
analysis, and column design principles. 

(R-22) - 6. 

(R-17) Rosenbrock & Storey - extensive coverage of practical mathemat- 
ical techniques for process dynamics. 

(S-8) Shreve - large number of industrial processes described in 

detail, with diagrams and non-technical descriptions, very little 
mathematics; book represents a lifetime of experience in the 
chemical industry; chemical engineering = unit processes 
(chemical changes) + unit operations (physical changes). 

(S-9) - 30. 

(S-3) Sawistowski & Smith - steady-state analysis and calculation. 
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(S-l) - 30. 

(T-l) Treybal - steady-state analysis and column design principles 
(B-12) Bird - mass transfer in terms of general conservation and 



diffusion equations; some transient analytical analysis. 

(H-5) Henley & Bieber - steady-state example, 3 tray column, a “ 2.5, 
benzene - toluene. 

(C-7) Chilton - cost analysis of distillation columns, y => cost 

$ installed per plate or foot, x = diameter squared, then the 
approximate equations for column cost ares 
Plate Columns - stainless - y = 6.4? 

Packed Columns - stainless - y = 9.82 ; 

distillation column maintenance 5~50 $/ft 3 bubble-plate, 

3-10 $/ft 3 sieve plate, 5-15 $/ft 3 for packed columns per 
year in I960 $. 

(M-2) - 43. 

(V— 3) Vilbrandt & Dryden - steady-state design of columns; "In 

distillation columns, the pivot point in design is the reflux 
ratio (B in this thesis), which can vary between minimum and 
total reflux. Higher reflux ratios require greater quantities 
of steam and cooling -water and a larger column diameter, but 
the column height requirements are lowered. The economic 
reflux ratio is usually 1.1 to 1.2 times the minimum for most 
cases," 

2. Extensive Bibliographies and Literature Surveys 

Literature Surveys from Industrial and Engineering Chemistry 
Process Control 

(w-16) (W-9) (W-15) (W-8) (¥-19) 
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(W-10) (W- 12 ) (R- 25 ) (W- 17 ) (F- 6 ) 
(W-ll) (W-13) (W-18) (G-10) (W-23) 



Distillation 

(B-17) (B-18) (B-19) (P- 8 ) (G-9) (F-7) 

(R-12)-I4ij (R- 8 )r?l. ; (H-7)-12.j (R-l 8 )-l 6 ,; (H-8)-10.; 

(R-17)-l.s (W-26)-4l . 5 (Z-3)-19. 

(G-14) Geddes - gives outline of developments that contribute to 

present knowledge of fractionator design in historical sequence; 
gives comments on present status and future problems; "If 
long-term funds were available for basic research on bubble 
plates, a substantial part of these should be assigned to 
scientific study of the fluid dynamics of plates." 

References of Historical Interest 

3. General Distillation 

(L-19) Lewis - very early (1909) presentation of basic principles; 
uses equilibrium and phase diagrams. 

(M-IO) McCabe & Thiele - the most referenced paper in the literature; 

original (L(25)McCabe-Thiele graphical diagram of steady-state 
column behavior; all previous methods (Sorel's was the first 
in 1899 ) analytical; algorithm for stepping off tray concen- 
trations on the equilibrium diagram. 

(L-25) Lewis - early ( 1922 ) presentation of stagewise steady-state 
calculations. 

(M-9) Murphree - explains McCabe-Ihiele Diagram in equation form. 

(V-2) Van Winkle - brief history of distillation - first recorded 

distillation in Egypt 50 BC, expect earliest unrecorded 
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distillation 2000 BC, fresh water distilled from sea water 
300 AD, beverage alcohol process first industrial distilla- 
tion process during 11-1 4th centuries, first books on 
distillation l 6 th century, stills were differential batch 
type with little reflux up until 19 th century, 19 th century 
began using steam (1800), bubble caps (1822), continuous 
still (1830), late 19 th century first recorded mathematical 
discussions of distillation by Sorel (1899) and Hausbrand 
(1893) » then 20th century began mathematics and improved 
distillation (L-19), (M-10), etc. 

(U-l) Underwood - the best single presentation of the history of 

distillation in the literature, many pictures and diagrams of 
ancient methods. 

(R-29) Rodebush - early (1922) plate-by-plate graphical technique; 
preceded McCabe-Thiele analysis. 

(E-3) Egloff - review and diagrams of 15th-l6th century distillation 
methods and apparatus. 

(C-14) Cope - early (1932) graphical method using McCabe-Thiele 
stepping type design on the lower part of the equilibrium 
diagram. 

(P-12) Peters - plate-by-plate calculations on the equilibrium diagram 
(1923) shortly after McCabe-Thiele. 

(T-6) Thiele & Geddes - referenced many times in the literature, 

graphical and mathematical steady-state analysis and calcula- 
tion method. 

(G-14) - 2. 
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(i-l) IBM 705 Program - (1959) steady-state computer program to 

perform Thiele-Geddes calculations; one of the first industrial 
programs. 

(A-17) Acrivos & Amundson - presents history, basic principles of 
distillation, steady-state analysis using numerical methods 
and eigenvalues. 

4. Dynamic or Transient Analysis 

(M-8) Murphree - (1925) first presentation in the literature of a 
discrete-plate dynamic equation; derives .holdup equation, 1 
ordinary differential equation, solves as an exponential. 

(B-29) Berg & James - (1948) early column transient behavior consid- 
erations; mainly concerned with startup problem, early presenta- 
tion of continuous-spatial model with partial differential 
equations, boundary conditions, and solutions, solves using 
linear equilibrium f(u) = mu + b, . experimentally verified results. 

(L-2) Lapidus & Amundson - (1950) early transient analysis; general- 
ized the results of (M-l) for countercurrent absorption, showed 
how outlet concentrations could be predicted from the time 
course of the two inlet compositions, assumes f(u) - mu + b 
equilibrium, uses Laplace transforms, poles & zeroes, linearized 
CSE, difference equations, entirely mathematical; calculates 
several transient responses. 

(B-7) Bartky & Dempster - (1948) early solution to the transient 

system of plate equations, analogous to the classic start-up 
problem except that top reservoir hats some holdup as individual 
stages, solution only approximate since compositions derived 
using a « 1.0. 
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(R-2) Rose & Johnson - (1953) early development of binary system model} 
relative volatilities, total flow rates, and holdups indepen- 
dent of time, discrete-plate equations solved by Euler predictor 
method, 

(R-12)-14.; (J-l)-21 . 5 (F-5)-19. 

5. Steady State Analysis and McCabe-Thiele Diagrams 

(E-4) Eckhart & Rose - steady-state prediction, discrete-plate 
equations, linear equilibrium f(u) = mu + b. 

(H- 26 ) Hartland - analytical, steady-state comparison, boundary 
conditions, 

(C-9) Cichelli - steady-state analysis, relates number of plates 
and reflux ratio to the sharpness of separation in binary 
batch distillation, graphs and equations for operation at any 
desired separation are given, many curves, sharpness of separa- 
tion defined in terms of the "pole height," Rayleigh equation 
used, 

(M-18) Mills - List of steady-state computer programs for equilibrium, 
enthalpy, etc, , calculations. 

(P- 8 ) Prausnitz - steady-state computer subroutines for column cal- 
culations for bubble temperature, dew temperature, bubble 
pressure, and dew pressure, 

(S-15) Surowiec - ideal cascade requires twice the minimum number of 
stages, 

(S-20) Strand - good development of discrete-plate steady-state 
equations. 

(Z-2) Zuiderweg - general steady-state comparison of device charac- 
teristics and plate efficiencies. 
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(F-3) Friday & Smith - discusses mathematically the formulation of 

a solution method for the equilibrium stage steady-state model; 
develops procedure for solving the concentration matrix equa- 
tions which avoids truncation error build up, does not require 
mesh points, works equally well for any number of feeds and 
side streams, and handles nondistributed components, 

(F-4) Furzer - nonuniform vapor distribution causes reduced efficiency, 
investigates plug flow model and perfectly mixed model, shows 
that maximum reduction in efficiency is halfway between these 
two models, 

(H-30) Himmelblau - general, recent process mathematical modeling 
text, block diagrams, frequency analysis, matrix analysis, 

(j-7) Jenson & Jeffreys - steady-state distillation analysis and 

mathematics, numerical methods for solving ordinary and partial 
differential equations, matrix methods, orthogonal function 
theory; primarily a mathematics text. 

(L-12) Lowenstein - uses normal (Gaussian) probability distribution 
scaled graph paper to greatly increase the accuracy of the 
McCabe-Thiele diagram at the ends, i.e. at very high and very 
low separations, 

(A-IO) Amundson & Pontinen - discrete-plate numerical steady state 

calculations, uses cubic equilibrium relationship, quadratic 
expression for enthalpy, temperature distributions, 4-5 
iterations on 15 plate column; results cited in (M-15). 

(B-26) Barker - efficiency, design, experimental analysis of bubble- 
cap trays. 
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(S-12) Sargent - steady-state multicomponent column numerical technique 
and digital computer design; considers individual plate 
efficiencies and equilibrium curve; no dynamics. 

(D-2) CEP Reprints - collection of steady-state distillation 

papers; articles covering heat and mass transfer, vapor- 
liquid equilibria, packed columns, tray column steady-state 
performance. 

(R-30) Rose - experimental evaluation of column steady-state. 

(S-24) Surowiec - steady-state column design using McCabe-Thiele 
diagrams and discrete-plate equations. 

(S-23) Stanislas - general steady-state characteristics 5 plate 
efficiencies, equilibrium, etc. 

(S-29) Sujata - plate-by-plate, steady-state calculations. 

(F-9) Floyd - locating feed trays for lowest cost operation, 
steady-state. 

(E-l) Edmister - true boiling point (TBP) defines distillation 

curve, distillation curve considered continuous and steady- 
state solved by graphical integration technique. 

(T-7) Treybal - graphical technique for finding plate efficiencies 
based upon McCabe-Thiele diagram, operating lines and equil- 
ibrium curve. 

(S-l6) Srygley - optimum steady-state design in the sense of minimum 
number of plates for desired product purity; uses Thiele- 
Geddes Q-method of convergence, sequential search for optimum, 
(B-l)-l.; (M-10)-3.;(H-10)-1.; (H-l8)-l. ; (M-8)-4. ; (M-9)-3.; 

(B-10)-30.; (C-l)-l. ; (V-2)-l.; (G-3)-l.| (H-5)-l. J (H-2)-l. ; 
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(A-17)-3.; (C-14)-3. ! (0-2)-l.j (T-6)-3 . j (R-5)-l6. ; (R-ZL)-l., 

(T-l)-l.j (S-3)-l.» (H-2l)-19. ; (l-l)-3.; (L-25)-3.» (A-19)-l . i 
( P- 12) -3. ; (R-l)-l 6 .j (R-29)-3.s (S-l)-30.; (S-7)-20. » (R-l 8 )-l 6 . 

6. Structural Design 

(L-13) Lowenstein - design of plate sizes using a nomograph and plate 
size "slide rule.” 

(R-22) Rose - minimum cost structural design considerations including 
wind loads, dead weight stresses, longitudinal stresses, 
thickness formulas, numerical examples, drawings of industrial 
columns . 

( J— 3) Jones & Van Winkle - experiments! analysis of 3 inch perforated 

plate column to determine plate thickness effects on column 
properties. 

(M-ll) Manning - structural screens introduced to provide better 
mixing and increase column efficiency. 

(L-12)-5. j (D-11)-16. 

7. Economics and Operations Analysis . 

(M-13) Mitten - economic optimization of distillation by dynamic 
programming. 

(B-15)-I.i (S- 8 )-l . 5 (F-9)-5.f (C-7)-l . i (V-3)-l. 

8. Thermodynamics 

(H-9) Hougen - Covers, very extensively, thermodynamics applicable 
to distillation. 

(R-23)-l. 

9. Hydrodynamics or Fluid Mechanics 
(F-4)-5. 5 (V-2)-l . t (B-12)-l.s (R-23)-l. 
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(L-26) Levlch - analytical solutions to convection and diffusion 
equations including chemical aspects. 

(H-14) Holm - concludes that vapor flow effects can sometimes cause 
Murphree efficiencies to be greater than unity. 

(P-5) -19. i (T-7)-5. 

(B-4) Bernard - considers sieve tray mixing, foam density, flow 
properties. 

(S-19) Sakata - time for mixing, analytical equipment, mixing pools, 
plug models, and tray efficiencies. 

10. Chemistry 

(B-l6) Black - simplified approach to phase equilibria, large number 
of phase equilibria for various compounds presented. 

(H-8) Hala - extensive list of chemical compounds referenced to 

entries in a large bibliography; many phase diagrams- presented, 
(B-l)-l.s (H-10)-l.; (F-9)-5.J (T-6)-3. 

(H-17) Harper & Moore - experimental paper showing a small still for 
measuring vapor - liquid equilibria lines, several lines given 
as examples. 

(H-20) Howard - concludes that any unsteady-state distillation cal- 
culations should include plate, condenses, and reboiler 
holdups in order to include enough degrees of freedom. 

11. Philosophy 

(A-13) Aris - oftentimes engineer’s rules of thumb are the only tools 
necessary to solve problems adequately, gives examples. 

14.2 DISTILLATION COLUMN DYNAMICS 

12. Textbooks 



- 52 - 



(M-l) Marshall & Pigford - beginning of quantitative analysis of 
unsteady-state operation; many simplifying assumptions used 
to obtain analytical solutions using Laplace transforms, 
assumed equilibrium relationship f(u) = mu + b with m, b 
depending only upon the identity of a component, assumed total 
flow rates and holdups independent of plate number and time; 
considers startup problem, 

(H-7) Holland - uses numerical methods to solve process differen- 
tial equations; uses 0-method for distillation discrete- 
plate equations, example step response of 3-tray column cal- 
culated; good literature survey; book emphasizes numerical 
methods and matrix methods, 

(G-3)-l. 

(F- 16 ) Franks - modeling of chemical processes for the purpose of 
control, plate equations, partial differential equations, 

13, Theses 

(M-15) Mohr - considers each of two sections of a column independently, 
determines step response for each input stream using IBM 705 
and discrete-plate equation, each of the response curves is 
then fitted with a two-time constant exponential expression, 
these expressions manipulated into general transfer function 
using Laplace transform variable p, time constants and gain 
factors of the column are then correlated with steady-state 
parameters from McCabe-Thiele diagram; analysis based on linear 
equilibrium f(u) = mu + b; assumes binary mixtures, constant cc» 
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constant fluid rates, total condensation, negligible holdup, 
constant liquid holdup each plate, reboiler holdup equivalent 
to that on each plate, q = 1, perfect mixing, unity plate 
efficiencies; found that major time constant was independent 
of the type of disturbance, gains and secondary time constants 
vary with disturbance type and with initial steady-state, 
predicted responses of large columns more sensitive to changes 
in values of selection parameters than small columns; column 
time constants increase with H/L, degree of separation, degree 
of nonlinearity of equilibrium curve, number of plates, 

(R-6) Romagnoli - hybrid simulation of discrete-plate equations; 

butadiene distillation plant, 2 towers of 49 trays each of 
bubble type, constant pressure, constant holdups, water 
ignored, vapor holdup ignored, constant Q, = 0.7, Francis- 
weir formula, micro assembly language and PDP-1 used; con- 
clusion was that hybrid setup worked faster than purely digital 
computation, 

(B-32) Brosilow & Tanner - develops and analyzes methods for computing 
and optimizing countercurrent staged processes using distil- 
lation as an example; formulates models in terms of both 
discrete and continuous spatial equations; solves continuous 
models using economic cost criterion using gradient and 
Lagrangian non-linear programming methods; found equations 
much easier to solve when models were not restricted to 
integral values of stages, allowing more general optimization 
methods to be used such as Lees’ method; two-point boundary 
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value problem solved by Imbedding it into an initial value 
problem in time; found that when Lees' algorithm was modified 
to generate tridiagonal matrices the rate of convergence 
increased substantially. 

(G-5) Gaydos - develops generalized digital computer program and 
model for a single bubble-cap plate; model includes multi- 
components, hydraulics, non-ideal vapor-liquid equilibrium, 
heat transfer to and from each stage, provisions for feed and 
take off streams; computation time limitations were encountered 
in the simulation. 

14. Reviews, Bibliographies, and Literature Surveys 

(A-3) Archer & Rothfus - presents survey of the 1955-1 9&0 dynamic 
behavior literature; discrete-plate equations developed, 
startup and transition between steady states, batch operation 
of plate and packed columns, and process control are all 
discussed in terms of their respective literature; frequency 
analysis. 

(R-12) Rosenbrock - surveys the history and present developments of 
discrete and continuous distillation and heat exchanges 
models; best presentation of this material in the literature, 

(H-?)-12. ; (R-8)-21. ; (Z-3)-19.f (L-3)-15.» (T-2)-l6. 

(S-33) Lehigh Symposium - discrete-plate models, transient response, 
and control of distillation columns; extensive bibliography of 
149 references on column steady-state, dynamic, and control 
aspects; frequency response analysis. 

Additional surveys, see 2. 
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Dynamic Models or Solutions 



Discrete Plate Equations 

15. Frequency Analysis or Laplace Transform Solution 
(A-3)-l4. ; (R-12)-l4. ; (H-7)-12. ; (R-8)-21. 

(Z-3)-19.i (S-33)-14.> (T-2)-l6. 

(L-3) Lamb, Pigford & Rippin - oscillations in tray compositions 

resulting from input oscillations in either feed composition 
or reflux flow are calculated using analog computer for 16- 
tray column; found frequency response at low frequencies like 
simple first-order process and at high frequencies having 
interference patterns and large phase lags; equations linearized 
about steady-state, E = 1, frequencies from 0.001 to 1 radian/ 
tray holdup time, 5-tray column frequency response also obtained. 

Transient or Time Analysis 

Numerical Solution 

16. Digital Computer 

(D-4) Distefano, May, & Huckaba - (1967), discrete-plate dynamic 

model solved for a sequence of upsets in which the next step 
occurs before the transients of the previous one have died out; 
solves large system of equations by Adams-Moulton-Shell (AMOS) 
finite difference technique on IBM 709, computation time was 
12 min. , expect hybrid steup would speed this up, 12-plate 
column, spacing 1 ft, 10" diameter, samples of every 3rd 
plate, 7 runs for different types of steps and pulses; 
oriented toward prediction for feedforward control; transient 
times 40-80 min. , experimental data agreed with computer 
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calculations to within 5 %t error blamed, on truncation in 
predictor-corrector methods, 

(B-32)-13. ; (R-4)-4. ; (R-20)-31. 

(D-8) Distefano - numerically solves discrete-plate equations by a 
large number of different methods and then compares them, 

(D-IO) presents stability aspects, 

(H-3) Huckaba - numerical solution by IBM 650 of transient response 
to binary 12-plate column, experimentally confirmed; uses 
nonlinear discrete-plate equations, inputs are step changes in 
heat input to the reboiler, feed composition, and simultaneous 
changes in feed composition and reflux ratio; computation time 
5 min, using fourth-order Runge-Kutta for starting and fifth- 
order modified Adams for continuing. 

(L-5) Luyben - uses set of linear perturbation type differential 
equations for characterizing dynamic behavior; experimental 
results verified for acetone -benzene system, 1,8 £ a < 2.2; types 
of disturbances were changes in feed composition, feed rate, 
top tray reflux rate, bottom tray vapor rate, equations 
solved by analog computer, 

(D-15) Duffin & Gamer - defines model for multicomponent distillation 
which includes secondary effects due to column hydraulics, 
holdup, and delay effects in boundary conditions; errors in 
generated system response at low values of elapsed time 
usually due to inadequate models; used general discrete- 
plate model and numerical integration scheme, conclude that 
model is valid. 
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(R-l6) Rosenbrock - surveys available solution methods for the discrete- 
plate equations, decides to use digital computer for speed and 
generality; describes numerical method and program, computa- 
tion time 5 min, for 5 plates up to 100 min. for 300 plates. 

(R-9) Rosenbrock - develops discrete plate equations, discusses 
possible methods of solution, almost identical to (R-l6), 

(R-10) Rosenbrock - discusses relative advantages of two computer 
programs for solving the equations in (R-9) » first program 
has f(u) fed in as table of 101 values, linear interpolation, 
forward integration used, first program finally rejected due 
to limitations; first program required solution of large 
system of simultaneous equations, second program developed to 
eliminate this by solving step-by-step by equating slopes; 
routine included to evaluate df(u)/du at discrete-points. 

(R-ll) Rosenbrock - discusses the accuracy of computer programs in 
(R-10), extends application to multicomponent systems; used 
equilibrium curve as f(u) = u + 0.02, E = 1.0; concludes that 
the most promising method for calculating transient response 
is digital computer. 

(A-4)-40, ; (R-17)-l. ; (S-33)-l^. S (S-l6)-5. 5 (S-20)-5. 

(M-3) Mah, Michaelson, & Sargent - dynamic behaviour of multistage 
systems described by large sets of non-linear first-order 
differential equations; discrete-plate equations then linearized 
but shown to be inconsistent due to linearization, step-by- 
step procedure using exponential function is proposed, 
numerical integration technique, reviews all standard numerical 
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(T-2) 



(W-l) 



(P-2) 



procedures, uses tridiagonal matrix, finds eigenvalues, shows 
that use of linear equilibrium f(u) = Ku also leads to physical 
inconsistency. 

Tetlow, Groves, & Holland - develops a generalized model 
which accounts for the effects of channeling, transfer lag, 
mixing, and mass transfer in unsteady-state, multicomponent, 
discrete-plate distillation; model tested on a large number 
of numerical examples; considers plug flow and perfect mixer 
holdups, large tri diagonal matrix results, 0 -method of con- 
vergence used; use of the generalized model in the analysis' 
of control systems is discussed. 

Wood & Armstrong - derives a linearized model of the discrete- 
plate equations using f(u) = mu + b equilibrium; solves for 
step response of feed composition; comparison with experimental 
results shows that the model is only valid for moderate values 
of time after the step and cannot be used as the column 
approaches final steady state. 

Peiser & Grover - presents a model similar to (H— 3) including 
the effects of heat and mass transfer and tray hydraulics; 
predicted that unsteady-state prediction can be used to solve 
several significant problems in multicomponent distillation 
which are not evident from steady-state analysis; numerical 
computations carried out on a digital computer simulating 
column open and closed loop control. 



- 59 - 



(D-9) Davison - solves large systems of x = Ax + Bu equations; 

calculates the poles and zeroes of the system, then solve 
for a few of the more significant variables in terms of poles 
and zeroes. 

(L-18) Luyben - discrete-plate equations solved for transient 
response for use in feedforward control. 

(D-6) Davison - discrete-plate equations solved using matrix methods 
for the transient response of a column due to pressure varia- 
tions for use in control. 

(S-6) Sargent - discrete-plate equations solved numerically using 
matrix methods. 

(T-9) Thorogood - discrete-plate equations solved using Runge-Kutta 
methods. 

(Y-l) Yesberg & Johnson - demonstrates use of a resistance network 
analog to solve the absorber problem of (A-l) and (L-2), 
discrete-plate equations linearized, time derivative represented 
by a backward finite difference to produce a set of simul- 
taneous algebraic equations which are solved on an IBM 650 by 
matrix inversion, 

(G-4) Greenstadt - discrete-plate equations solved by Newton's 
method, 

(R-5) Rose, Sweeney & Schrodt - solves startup problem for ternary 
mixture using discrete-plate equations, Lewis-Matheson method 
used. 

(D-ll) DiLiddo & Walsh - pulse column considered as a series of 

stages, 3 plate ideal model, 9 plate real model, numerical 

solution on IBM 605. , 
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(R-l) Rose, Johnson & Williams - (1950) discrete-plate equations 
solved by Euler predictor method; model assumes relative 
volatilities, flow rates, and holdups independent of time, 

(R-19) Rose, Johnson & Williams - (1951) similar to (R-l), early 

papers showing pictures of IBM cards, plate equations solved 
by numerical methods, 

(L-27) Lowe - discrete stage equations solved on digital computer, no 
distillation, 

(W-4) Waggoner & Holland - discrete-plate equations solved using 

Simpson's rule, 3 point corrector to approximate the integrals 
in the component material balances, results presented in table 
form. 

(R-18) Rose & Williams - (1950) early plate-by-plate solution of 
discrete model, similar to (R-19) (R-l), shows wiring of 
computer control panels, 

17. Hybrid Computer 

(F-l) Franks - solves several countercurrent problems, several of 
these were for distillation. 

(F-13) Frank & Lapidus - discrete-plate equations solved by hybrid 
computer using one integrator for each plate, 

(R-6)-13. j 

(F-12) Frank & Lapidus - hybrid simulation particularly useful for 
nonlinear partial differential equations, uses repetitive 
analog operation, memory, and first-order lags. 
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18. Analog Computer 



(B-23) Bowman & Clark - uses linear equilibrium f(u) => ku which was 
found to be valid near the top of the column where volatile 
component concentration is small, discrete equations for 20 
and 30 plate columns wired up on analog computer; when column 
was switched from total reflux to finite reflux ratio, an 
almost instantaneous drop in overhead, composition occurred, 
from this point on the overhead was independent of the time 
on total reflux and dependent only on the stillpot composi- 
tion at the end of the total reflux period; this suggests that 
column can be divided into two independent sections, first 
would be period on total reflux, second would be behavior at 
finite reflux; linear equilibrium does not provide exact rep- 
resentation of real column, provides guide to the degree and 
direction of change within the column. 

(P-l) Pigford, Tepe & Garrahar. - solves unsteady-state equations for 
batch distillation of a binary mixture using a mechanical 
analog computer called the "differential analyzer" ; assumes 
constant relative volatilities, vapor & liquid flow rates. 

(R-3) Rose & Williams - demonstrated use of an analog computer to 

obtain transient response and design a controller for a 5-plate 
column, uses analog computer with Pads delay circuits; large 
number of problems solved to determine the best controller for 
maintaining the composition of the distillate constant under 
composition and thermal variations of the feed. 

(G-12) Grover & Peiser - analog computer solves plate equations for 
control, stability aspects considered. 



(L-5)-l6 
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19. Analytical Analysis 



(D-3) Davidson - uses mechanical analog and Rayleigh's method 

to determine eigenvalues for finite series representation of 
plate column transient behavior; plate-type stripping column 
fed at top, no bottoms takeoff, assumes linear equilibrium; 
equilibrium concentration on any plate proportioned to 
exp (-0T), where g depends upon a and number of plates N; 
Rayleigh method used to give approximate g; solves example 
model of Taylor Diffusion type, found that first term in 
series was most important. 

(G-2) Gilliland & Mohr - analytical analysis of discrete-plate 

equations using two-time constant exponential model of (M-15) ; 
digital computer solves for transient response, then the two 
time constants are determined from the results and used to 
develop transfer functions; by use of transfer functions, 
the responses of several columns to step changes in feed 
composition were predicted and compared with the responses 
calculated by numerically solving the discrete-plate equations. 
(A-3)-l4.; (M-8)-4.; (M-l)-12. ; (R-8)-21. ; (S-3l)-43. 

(R-7) Rosenbrock - funcamental "disturbance trapping" paper; char- 
acterizes the departures of a distillation column from its 
steady-state by a quantity D which measures the rates of change 
of composition on all the plates and increases whenever the 
column is disturbed from its steady-state; result is generalized 

D = 2 ?| d(H n u „)/ dt | 
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(R-32) 

(W-2) 

(F-16) 

(3-27) 

(C-2) 

(Z-3) 

(F-5) 



to multicomponent systems; disturbances become "trapped; in 
the column at points where equilibrium curve slope equals 
operating line slope, rates of change of composition do not 
decrease as would be expected, numerical computation of behavior 
in these cases is difficult; presents an electrical analog 
to the discrete-plate equations. 

Rosenbrock - discrete-plate equations, matrix methods, energy 
considerations . 

Wilkinson & Armstrong - considered response of a column at 
total reflux to a change in bottom vapor composition, used 
binary mixture, linear equilibrium curve; response predicted 
by analytical analysis was in good agreement with that observed 
experimentally. 

■12.; (Q-3)-l.S (H-20)-10. ; (3-7 ) -4. ; (B-12)-l. 

Balasubramanian - analytically solves one differential equation 
for a one-plate still, 

Cullinan - considers transient start-up problem using matrix 
methods for analytical solution, uses f(u) = mu + b. 

Zykov - analytical solution of discrete-plate equations for 
transient analysis of multicomponent column; good bibliography 
of Russian literature. 

Foss, Gerster & Pigford - assumptions of complete mixing or 
no mixing lead to inaccuracies in distillation, paper attempts 
to establish the nature and extent of mixing and to develop 
calculationai methods to account for its effect on plate 
efficiency; mixing experimentally determined by use of 

tracers and measurements of residence times; simplified calcu- 
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lation procedure presented which affords a rapid means of 
computing plate efficiency under all mixing conditions. 

(H-2l) Hassett - solves for transient behavior using equilibrium, 
steady-state, and McCabe-Thiele diagram. 

Continuous Spatial Equations 

20. Frequency Analysis or Laplace Transform Solution 

(J-6) Jafri, Glinski & Y7ood - continuous system transient response 

with control using time constants and transfer function analysis. 

(M-27) Majumdar - solves CSS using Laplace transforms and lineariza- 
tion, applies to Clusius column. 

(H-4) Hoerner & Shiesser - frequency and time responses using Laplace 
transforms, linearized CSE, linear equilibrium; gives model of 
Taylor Diffusion type. 

(S-7) Sellers & Augood - transients in a liquid hydrogen separator, 
uses Laplace transforms, exponential characterizations, rate 
of approach method. 

(W-5) Ward - uses Laplace transforms and frequency analysis to cal- 
culate the time behavior of dynamical systems. 

(D-l) Douglas & Eippin - uses linearized equations and sinusoidal 
inputs, considers system in terms of chemical oscillators, 

21. Transient or Time Analysis 

(j-l) Jackson & Pigford - (1956) digital computer solution of linear- 
ized CSS model for startup problem, plots of composition 
throughout column as a function of reduced time are presented; 
linearized CSE and linear equilibrium f (u) = Au gives equation 
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of Taylor Diffusion model type; took 3 hours of computation 



time on IBM 701 to solve for transient curves on several trays. 
(H-4)-20, ; (D-3)-19. ; (M-l)-12. ; (C-2)-19.; (L-26)-9. 

(K-2) Koppel - solves heat exchanger/chemical reactor equation of a 
form possibly similar to CSE; 

9u = £l + r(t)~j 3u + PjjL + br(t)"| u n 

at 

(O-l) Osborne - linear equilibrium f(u) = mu + b used, equations 

considered continuous in time and differenced in theoretical 
stage; concludes that the real cause of numerical instability 
problems is in theoretical stage direction differencing; 
numerical solution obtained with very large theoretical stage 
step, one theoretical stage, and very small time steps, this 
took care of instability problems; Fortran IV computer program 
for max. 6 components, max. 60 trays described, 

(R-8) Rosenbrock - calculates transient responses; describes labor 
needed to solve CSS models; describes control aspects and 
theory; presents good bibliography. 

(S-13) Stone & Brian - detailed numerical methods described for 
solving CSS type equations; solves Taylor Diffusion model 
equation as an example; CSE is a subcase of equation (2) 
which is a general form of convective transport equation; 

[D(x,t,u) £u] - d__ [V(x,t,u)f(u)] = £u 

ax a x 

several different types of numerical methods discussed and 
compared; analysis derived applies only to linear equations; 
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recent results indicate that the desirable features of 



solutions obtained by the new equations for linear problems 
are to a large degree found in solutions of nonlinear problems. 

(M-28) Montroll & Newell - exact solution of nonlinear differential. 

equations which describe the time dependent behavior of multi- 
stage cascade separating processes of two very similar non- 
linear species, Rayleigh separation law postulated for each 
stage? linearization of nonlinear second-order portion differ- 
ential equations is discussed; analytical expressions in terms 
of exponentials and eigenvalues are developed; equilibrium 
curve is approximated by f(u) = u + cu(l - u). 

(K-3) Kermode & Stevens - several nonlinear continuous models 
solved on an analog computer, 

(P-7) Powers - numerical, solution to Taylor Diffusion type partial 
differential equation with two-point boundary conditions and 
one initial condition, 

(R-14) Ruckenstein - analytical analysis of convective diffusion 
(Taylor Diffusion) type equations; extensive analysis of 
applicable transformations and linearization methods, 
(B-29)-4.; (G-3)-l. 

(T-5) Tsang - solves heat equation using eigenfunction expansion 
(orthogonal), evaluates 1st 3 values numerically; other 
analytical solutions presented also, Bessel functions, 

(H-25) Herron & Van Rosenberg - uses a mesh and centered difference 

method to numerically solve convective transport equations with 
two-point boundary conditions. 
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(j-2) Jury - solves heat equation using analog computer, uses 
memory, solves repetitively. 

(J-5) Jackson - numerical solution and optimization of partial 
differential equation models. 

(L-20) Lapidus - transient response and control of chemical reactors 
using continuous models. 

(B-20) Bedingfield & Drew - heat and mass transfer expressed by the 
same equations. 

(W-14) Woodle - analogy between distillation and heat transfer, some 
equations. 

(B-25) Brian - transient response using a continuous, Taylor 
Diffusion type model. 

(C-8) Crank - diffusion in different geometries with different 
boundary conditions, finite difference methods, diffusion 
and chemical reaction; use transformations for equations with 
variable diffusion coefficients. 

Experiments! Transient Behavior 

22, Frequency Response 

(H-15) Henley - frequency response techniques for experimental 
analysis of transient behavior. 

(H-23) Hutchinson & Shelton - frequency response techniques using 
correlation functions. 

(A -6) Armstrong & Wilkinson - carried out experimental work to verify 

theoretical computations and methods of (R-l6) and (R-9) 5 
studied behavior of 21 plate CH^ - CCl^ separator subject to 
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step changes in feed and reflux on column composition; long 
time agreement was better than short time after disturbance; 
results summarized by transfer functions having the form of 
pure time delay followed by a linear log, both time constants 
are functions of plate number. 

(H-24) Haagensen - experimental results derived from frequency 
response using matrix techniques. 

(W-28) Woods - experimentally determined controller settings for a 
continuous system operating dynamically. 

23. Time Response 

(H-3)-l6 . 5 (L-5)-l6.; (R-30)-5. 

(B-5) Baber & Gerster - determines experimental transient response 
of column to changes in liquid and vapor rate, demonstrate 
the applicability of discrete-plate equations for predicting 
measured response; responses to step inputs presented in 
tables and graphs; equations solved by analog computer; 
results confirmed validity of model; linear perturbation types 
of equations predict satisfactorily the transient behavior. 

(B-3) Baber - similar to (B-5) but earlier; experimental response 

of 5 “tray, 2 ft, bubble-cap column; tests made over a range of 
gas and liquid rates nearly up to the flooding point and for 
tower pressures up to 5 atmospheres; average difference between 
predicted and experimental results was 13^> indicating that 
simple perturbation equations are valid models, 

(R-15) Rademaker - tested an ethylene-ethane splitting column to 

provide data for checking a general theory of column dynamics; 
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data summarized, in graphs and. experimental accuracy discussed,. 

(D-12) Davies - "hidden transients" can have significant effects on 
tray efficiencies. 

24. Cyclic Distillation 

(A-14) Atkeson - one of the first papers (1957) to indicate that 

cycling or nonequilibrium operation can increase mass transfer 
rate greatly; attempts to explain physically. 

B2.3 DISTILLATION COLUMN CONTROL 

25. Textbooks 

(A-12)-40.; (G-3)-l. {(C-l)-l. 

(B-13) Buckley - very practical, industrially oriented process control 
text. 

(K-5) Koppel - matrix theory, optimal control, sampled-data control. 

(A-ll) Athans & Falb - mathematical theory of optimal control; many 
examples and problems; very clear presentation; extensive 
bibliography. 

26. Theses 

(B-2) Beecher - presents method of dynamic control system synthesis 

using calculus of variations; provides extensive Fortran program; 
example is 100 plate butadiene/butiene-2 column. 

(B-32)-13. 

(G-l) Gordon-Clark - uses matrix theory to adjust dynamic response 
of process to some required response before applying conven- 
tional control; uses linearized model of a 5-pla.te binary 
column, control and measure composition on each plate; main 
computation difficulty is finding eigenvalues of matrix; 
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Fortran program, form matrix, form into tridiagonal matrix, 
obtain complex roots of real polynomial by Bairstowe's method, 
iteration using quadratic factor, obtain roots of quadratic 
factors, return; major drawback is large number of variables 
required. 

27. Extensive Bibliographies and Literature Surveys 
(A-3)-l4 . 5 (R-8)-21.j (R-12)-l4.; (S-33)-l4. 

Conventional Control Systems 

28. Digital Control 

(H-29) Hanson, Duffin & Somerville - provides extensive Fortran 
computer programs for control. 

(B- 6 ) Buster - closed loop control applied to oil fractionating 
system, required large memory and cleyer programming to 
complete calculations in real time. 

29. Hybrid Control 

see 17. 

30. Analog Control and Instrumentation 

(H-16) Haines - large number of diagrams of possible column control 
configurations with explanations, 

(G-3)-l.j (B-13)-25.J (C-l)-l . 5 (G-12)-18. 

(B-22) Buckley - basic column control strategy, diagrams of control 
loops, linearized models. 

(B-IO) Bauer & Orr - uses McCabe-Thiele diagram to derive best 
operating lines for control. 

(B-31) Boyd - mostly nonmathematical and non-diagramatical discussion, 

(P-6) Pink - describes control using analog computer, no mathematics. 
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(s-1) 



(L-4) 



(C-12) 

(T-ll) 

(M-29) 

(P-14) 

(R-13) 

(S-9) 

(H-l) 



Shinsky - process controls based upon time response analysis; 
distillation columns hard to control because (l) many trayed 
towers slow to respond to control action, (2) separation 
requires many variables, (3) on-line analysis not always 
available, (4) distillation units are usually last in the chain 
of processing operations, and (5) factors affecting separation 
not readily interpreted in terms of control system requirements, 
Lupfer & Parsons - describes a control system designed to reduce 
the effects of changes in flow rate and feed composition on 
column operation, uses prediotive or feedforward control, 
dependent outputs of the process are controlled by measuring 
one or more inputs (which generally cannot be controlled), 
and then the controlled parameters are changed as required 
to achieve the desired output, 

Ceaglske - comparison of transient and frequency response 
methods for control of linear chemical process systems, 

Tivy - control discussion, no mathematics. 

Moczek - discusses effect of transient behavior and dead time 
on control, very little mathematics, 

Phillips - control descriptions, no mathematics. 

Rijnsdorp - discusses feedforward & feedback control with 
distillation column as an example. 

Strobel - describes in detail theory of optical and electrical 
measuring devices; great deal of theory presented, photometers, 
spectral analysis, wave theory, etc. 

Harriott - transient and frequency analysis used for column 



control 
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Control Systems Using Dynamic Models 



31, Digital Computer Control 

(A-12)-40. ; (D-4)-l6. 5 (L-l8)-l6. ; (D-6)-l6. 

(R-20) Rosenbrock - application of automatic control theory to 

chemical processes has not led to the same improvement of 
performance as it has in the control of mechanical and elec- 
trical systems, paper attempts to explain why? complexity of 
chemical processes, lack of suitable measuring equipments, 
remote objectives i.e, out of 1000 variables only 10 are of 
interest, modes not always easily separated; proposes matrix 
method in which important modes are picked out for control, 
(S-33)-l4. 

(C-ll) Cadman, Rothfus & Kermode - uses matrix methods and frequency 
response techniques for design of multicomponent feedforward 
control system, 

32, Hybrid Computer Control 

(D-13) Dahlin & Nelson - uses hybrid computer and matrix maximum 
principle for optimal control. 

33, Analog Computer Control 

(j-6)-20. ; (K-2)-21. ; (L-3)-15.J (R-3)-l8.; (H-15)-22. 

(L-14) Lupfer & Oglesby - elaborate analog controller described, 
instrumentation and feedback schemes, 

(L-6) Luyben & Gerster - studies effectiveness of feedforward 

control for 10 and 40 tray columns, performance of overhead 
and bottoms controller determined by analog simulation and by 
experimental tests on a 10 tray - 2 ft, column; concludes that 
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relatively simple feedforward controllers appear adequate for 
distillations for small input disturbance, a linear model can 
be used to determine controller transfer functions, 

(W-27) Williams & Harnett - uses frequency response analysis, plate 

equations, first order lags; describes various control schemes. 

34. Optimal Control 

(L-9) Lapidus - nonlinear optimal control, quadratic performance 
criteria, Riccati equation and solution, etc. 

(B-ll) Brosilow & Handley - optimal control of the overhead composi- 
tion of a distillation column, integral squared error criterion 
on disturbances; experimental analysis on 5 inch column with 
15 trays and 3 bubble caps per trayj control system behaved 
well in spite of model inaccuracies. 

(B-32)-13.; (A-11)-25.S (D-13)-32.s (K-5)-2 5.i (S-33)-l4.s 

(j-5)-21. 

35. Distributed or Modal Control 
(D-9)-l6. ; ( J— 6) — 33- t (S-3l)-43.s (F-l6)-12. 

(G-13) Gavalas - eigenvalue solutions for distributed parameter 
steady state. 

B2.4 MATHEMATICS AND COMPUTATION 

36. Ordinary Differential Equation Theory 

(H-ll) Hartman - $20.00 text, ordinary differential equations in 

all aspects from a pure mathematics standpoint, theorem-proof 
presentation. 

(B-14) Birkhoff & Rota - very good presentation of transformations and 
eigenfunction expansions for two-point boundary-value problems 

and Sturm-Liouville problems. 
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Ince - if one were ever restricted, to only one look on ordinary 
differential equations, this would have to be the one. 

(H-12) Hildebrand - very practical applied mathematics text, useful 
theory of infinite series expansions to solve DS's. 

3?. Partial Differential Equation Theory 

(G-7) Garabedian - theory of 1st and 2nd order PDE's, basics of 
integral equation theory, very few examples or problems. 

(F-ll) Forsythe & Wasow - practical numerical methods text, useful 
for solving PDE's. 

(B-2l) Berg & McGregor - very clear presentation of basic principles 
and solution techniques for PDE's. 

(0-1) -21. ; (F-12)-17. ; (F-l6)-12. 

(G-6) Gurel & Lapidus - extensive discussion of stability of ODE's 
and PDE's. 

(W-?) Webster - large number of practical example-problems solved. 

(H-28) Hildebrand - newest Hildebrand text, seems to be as practical 
and understandable as the previous texts; numerical methods 
for solving PDE's. 

38. Integral Equation Theory 

(T-3) Tricomi - best text on integral equation theory in literature, 
complete, readable, very few examples. 

(P-4) Petrovskii - good presentation of separable kernel theory and 
use of algebraic equations to approximate integral; several 
classifying examples and problems. 

(L-8) Lovitt - very clear presentation of basic theory with large 

number of easily worked and instructive examples and problems. 
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(H-13) Hildebrand - very explicit presentation of methods to convert 
from IE's to DB's with examples. 

(D-5) Davis - large number of methods and equations solved which 
are not found anywhere else in the literature. 

(M-6) Mikhlin - very readable text, terminology sometimes different 
from standard U.S., coveres about same material as (L-8). 

(S-ll) Smithies - rigorous text, presents theory of IE's in terms 
of Lebesgue integration and L 0 spaces; Riemann integration 
and spaces are subcases of L 2 spaces. 

(V-4) Volterra - original classic in theory of integral equations. 
(G-7)-37. ; (W-7)-37. 

(W-6) Whittaker & Watson - very clear, but abbreviated, presentation 
of IE theory; extensive presentation of infinite series 
expansions. 

39. Mathematical Transformations 

(T-4) 'Tranter - basic theory and application of different types of 
integral transforms. 

(Z-l) Zemanian - recent text, extensive theory of transformations 
and integral transforms. 

(S-3l)-43.; (m- 16)-44. ; (R-l4)-21.; (C-8)-21. 

40. Matrix Mathematics 

(A-4) Acrivos & Amundson - transient solution to column discrete 

equations by matrix methods, e^ solution; this paper refer- 
enced very often in the literature; paper originally brought 
the subject of matrices to the attention of chemical engineers, 
wide variety of chemical engineering problems solved in this 
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paper by matrix methods including the absorber problems of 
(A-l) and (L-2). 

(A-12) Amundson - good presentation on eigenfunction expansions; 
matrix theory. 

(h-13)-38.; (J-7)-5. 

41. Numerical Solution Techniques 

(B-30) Berry & DePrima - develops iterative procedure for the 

determination of the eigenvalues and eigenfunctions associated 
with the solution of Sturm-Liouville problems in a finite 
interval; presents and discusses convergency of an iterative 
scheme different from "sweeping" or Rayleigh-Ritz; presents 
numerical example. 

(R-24) Ralston - very practical and widely referenced text, many of 
the programs in the Scientific Subroutine Package (I— 3) use 
methods of this text. 

(D-8)-l6.; (H-28)-37. ; (F-ll)-27. 5 (B-32)-13. 

(F-2) Fox - somewhat dated but clear presentation, verified with 
many numerical examples, of methods for finding boundary- 
value solutions and eigenvalues. 

(K-l) Kenneth & McGill - gradient methods, general numerical methods, 
some theorems on existence and uniqueness of boundary-value 
solutions. 

(L-10) Liu - numerical solution by finite-difference method; shown 
to be explicit, stable, more accurate than Crank -Nicholson; 
solves heat equation and several nonlinear examples. 
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(M-22) McGinnis - solves Taylor Diffusion type equations using 

Runge-Kutta techniques; describes numerical methods for BV 
problems. 

(S-13)-21.; (C-8)-21. ; (A-10)-5. I (P-8)-5. 

(H-19) Hamming - very practical and extensive numerical methods text. 

(M-4) McCracken & Dorn - application of numerical methods, flow 
charting, basic Fortran, programming principles, ODE's and 
PDE’s numerical, methods, etc. 

(R-4) Rose, Johnson & Williams - used both analog and digital 

computer to solve plate equations for a 7 plate binary column; 
found that the time required for the column to change from one 
steady-state to another after an abrupt change in feed composi- 
tion is a strong function of the magnitude of change, a, reflux 
ratio, N, feed tray location; numerical methods used described 
briefly. 

(R-29)-3.S (P-12)-3. ; (H-25)-21.; (J-7)-5. J (A-17)-3. 

(L-ll) Lee - invariant imbedding approach; classical methods use 

quasi-linearization; invariant imbedding considers a family 
of problems from zero to the duration of the original problem, 
by imbedding, solve for the missing conditions for 2-point 
BV problem; solves Taylor Diffusion type equation as an example. 
(D-9)-i6. 5 (m-2)-43.i (R-l7)-l. ; (J-5)-21. 

(N-l) Naylor - numerical solution techniques. 

(W-26) White - method for numerically calculating eigenvalues and 
eigenvectors of large dimension matrices; good bibliography. 
(D-15)-l6.; (S-12)-5. 
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42. Boundary Value Problems 



(B-8) Boyce & DiPrima - basic text on practical solution techniques 
and theory of boundary and eigenvalue problems. 

(B-9) Beltrami & Wohlers - very theoretical presentation; existence 
and uniqueness theorems for boundary value problems. 

(V-l) Villadsen & Stewart - new collocation methods given for 

solving symmetrical boundary-value problems using orthogonality 
conditions to select collocation points; accuracy is shown to 
be comparable to least squares or variational methods, 
calculations are much simpler; applications given to one- 
dimensional eigenvalue problems and to parabolic and elliptic 
PDE's; collocation methods are special techniques for solving 
integral equations numerically. 

(L-ll)-4l. ; (H-ll)-36.; (K-l)-41. ; (B-l4)-36. ; (M-22)-4l. ; 

(B-32)-13. ; (F-2)-4l. ; (H-12)-36. ; (H-25)-21. 

43. Eigen-Values , Vectors, and Functions 

(S-3l) Singer - state variable transformations and matrix methods 

to select significant modes and eigenvalues of multivariable 
systems. 

(M-2) Mickiey, Sherwood & Reed - numerical solutions of PDE's 
using finite differences applied to stagewise processes; 
orthogonal functions and infinite series solutions of PDE's. 
(B-30)-4l. ; (A-12)-40.; (G-13)-35 . i (H-ll)-36. ; 

(H-13)-38.j (l-2) - 36 . ; (B-l4)-36. ; (D-3)-19. 5 

(S-13)-21.; (V-l)-42, ; (W-26)-4l . 5 (B-5)-22.; 

(H-12)-36.; (J-7)-5. 
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44. Special Functions 



(A-18) Abramowitz & Stegun - complete and extensive analytical , 

numerical, and graphical presentations of special functions. 

(L-23) Lebedev - practical user's text showing mathematical properties 
of special functions and examples of usage. 

(M-I 6 ) Magnus, Oberhettinger & Soni - states and proves many of 
the mathematical properties of special functions. 

(R-28) Rainville - similar to (L-23). 

45. Computation and Computer Programming 

(0-4) Organick - basic principles of Fortran IV in textbook form 
with large number of examples. 

(l-3) Scientific Subroutine Package - set of over 250 subroutines 
for performing standard numerical manipulations such as 
matrix inversion, integration, differentiation, expansion in 
functions, least squares curve fitting, roots of polynomials, 
etc. ; (l-3) is in Fortran IV but recently ( 1969 ) a reduced 
package in PL/l has been released. 

(M-4)-4l. ; (P-8)-5 . 5 (K-29)-28. ; (M-l8)-5. ; (B-15)-l. 
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SECTION 2 



MODELS OF BINARY DISTILLATION COLUMNS (m) 
Ml THE DISCRETE-PLATE EQUATIONS (DPE) 

M2 THE CONTINUOUS-SPATIAL EQUATION (CSE) 
M3 SOLUTION TECHNIQUES FOR THE CSE 
M4 LINEAR APPROXIMATIONS TO THE CSE 



"I HAVE HARDLY EVER 
OF REASONING.” - PLATO 
"PLATO WAS A FOOL!" 



KNOWN A MATHEMATICIAN WHO WAS CAPABLE 



- JACKIE GLEASON 



THIS SECTION PRESENTS DISCRETE AND CONTINUOUS MODELS FOR THE 
COMPOSITION BEHAVIOR OF A BINARY PLATE DISTILLATION COLUMN. THE 
CONTINUOUS-SPATIAL EQUATION (CSE) IS DEVELOPED FROM THE DISCRETE- 
PLATE EQUATIONS (DPE) BY CONSIDERING THE PLATE NUMBER TO BE A 
CONTINUOUS VARIABLE. SOLUTION METHODS FOR THE CSE ARE PRESENTED, 
INCLUDING TRANSFORMATIONS, GENERAL SOLUTION TECHNIQUES, NONLINEAR 
APPROXIMATIONS, AND LINEAR APPROXIMATIONS. AS A FINAL RESULT THE 
LINEAR POLYNOMIAL-COEFFICIENT MODEL (LPCM), WHICH IS THE SUBJECT OF 
SECTION 3, IS DEVELOPED AS A LINEAR APPROXIMATION TO THE CSE. 
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CHAPTER Ml 



THE DISCRETE-PLATE EQUATIONS (DPB^ 

The development of the discrete dynamic model of a binary plate 
distillation column is the subject of this chapter. The development 
begins by defining the symbols to be used in the derivation. Then 
the individual plate equations will be developed by considering com- 
ponent mass balances on a typical plate. Finally several comments are 
made referencing the literature pertaining to the discrete-plate 
equations. The developments of this chapter are a simplified version 
of those presented in reference (G-3). 

The symbols to be used in the derivation of the individual plate 
equations are presented in Table Ml.l. These symbols are to be 
applied to the characteristics of the idealized, typical plate of the 
lower or stripping section of the column as shown in Figure Ml.l, 

The derivation of the equation for the rectification or upper section 
are identical except that the upper flow constants must be used. The 
author has attempted to keep these symbols consistent with those used 
in Chapters II and 12 and especially in Table 12,1, 

Ihe dynamic model is a system of equations developed for each 
plate by utilizing the basic principle of conservation of mass stated 
in Equation Ml.l. The four possible mass balances which can be applied 

Accumulation = Inflow - Outflow Ml.l 

to the plate of Figure Ml.l are listed below. Ihree of these mass 
balances, numbers 2, 3» and 4 are satisfied by the assumptions stated 
in Table Ml.l for the individual variables. These assumptions restrict 
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u - Concentration of the lighter component in n-th plate liquid 
(lbm mole lighter component/lbm liquid) 

f(u n ) - Concentration of the lighter component in the n-th plate vapor 
(lbm mole lighter component/lbm vapor) 

a - Relative volatility (dimensionless); assumed constant 

F - Feed rate (lbm liquid/hour) ; assumed constant 

D - Distillate rate (lbm liquid/hour) ; assumed constant 

W - Bottoms withdrawal rate (lbm liquid/hour) ; assumed constant 

H - Liquid holdup on the plate (lbm liquid) ; assumed constant 

h - Vapor holdup above the plate (lbm vapor) ; assumed constant 

L - Liquid rate in the rectification section (upper section) 
u (lbm liquid/hour); assumed constant 

L, - Liquid rate in the stripping section (lower section) 

(lbm liquid/hour) ; assumed constant 

V - Vapor rate in the column (lbm vapor/hour); assumed constant 

- Upper reflux ratio, L u /V, (lbm liquid/lbm vapor) 

B, - Lower reflux ratio, L^/V, (lbm liquid/lbm vapor) 

k - Feed plate index (integer) 

n - Internal plate index (integer), 1 < n < N 

N - Total number of plates in the column (integer) 

E - Murphree plate efficiency (dimensionless), defined by 
Equation 12,1, E = 1 assumed for all plates 

q - Portion of the feed which adds to the lower liquid rate; 

L = L + qF ; further developments use q = 1, saturated 
liquid U feed 

Table Ml.l - LIST OF SYMBOLS USED IN DISTILLATION 
COLUMN DYNAMIC MODELS 
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1 



. Liquid phase mass balance of the lighter component, 

2. Liquid phase overall mass balance, 

3. Vapor phase mass balance of the lighter component. 

4. Vapor phase overall mass balance, 

the overall validity of the model greatly but are utilized to simplify 
the developments which follow in later chapters. It is eventually 
expected that the generality of the Linear Polynomial -Coefficient Model 
(LPCM) presented in Section 3 can be utilized to generate solutions 
which are nearly as accurate in representing the essential nature of 
the column transient behavior as would be a completely general model 
using four equations per plate. 



i 




Figure Ml.l - TYPICAL RECTIFICATION PLATE WITH SYMBOLS 

33P 



The mass balance of the lighter component in the liquid phase 
gives equation Ml, 2 for the upper section and Ml. 3 for the lower 
section. The steady-state portions of these equations are seen to 
be equivalent to the discrete-plate steady-state equations in 
Figure 12, 3« 

H du n 

d — - V [ f K-l> - f (“„)■) + L u [“n+l-"„3 m.2 

H du n 

dT = V WVl) “ f(u n)l + L lC u n + l- u n l ^*3 

The overall binary distillation column dynamic model utilizing 
the discrete -plate equations can now be presented as Figure Ml, 2, 
where the effects of the condenser and reboiler have not been included. 
It can be seen that this model consists of N-nonlinear ordinary differ- 
ential equations with boundary conditions. This particular model is a 
highly simplified version of the more general models (using all 4 
equations with varying holdups and flow rates) usually considered in 
the literature (See B2»2 - Discrete Plate Equations), The solution 
techniques to be applied to the continuous-spatial equation developed 
in later chapters can also be applied to more general models; but for 
the purposes of presenring the LPCM and its solution techniques, the 
model of Figure Ml. 2 is sufficient. When accurate quantitative behavior 
is desired, the discrete-plate model is almost always the model 
employed in the control system applications described in the literature 
(B2.3 - Digital Computer Control). 

The reader is referred to the very extensive literature in 
Chapters B1 and B2 for further information pertaining to discrete- 
plate models and solutions. This chapter has presented a very brief 
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development of the discrete-plate model (DPM) for the purpose of 
utilizing it in the development of the continuous-spatial model (CSE) 
which is subsequently to be linearized to the linear polynomial- 
coefficient model (LPCM) for which an analytical solution procedure 
is developed. The next chapter presents the development of a contin- 
uous-spatial model (CSE) from this discrete-plate model. 



Model Equations 

Rectification Section k + 1 < n < N 

air = v C f K-i) ' f(u n)l + L u IX+i- u n l 



Feed Tray 
dt 



n * k 



- V C f ( u k-l) - + L u u k+1 - L l u k + Fu f 



Stripping Section 



1 < n < k - 1 



H du 



dT = v L f K-i) - + L i LViAl 



Equilibrium 

f(uJ - 



au. 



1 + (a-i) 

End Conditions 



f(u ) = u , 
n d 



U 1 “ S 



Top n 
Bottom n 



N 

1 



Feed and Flow Conditions 

= L u + qF or B 1 = B u + qF/V 

V = L + D 
u 

W = L 1 - V 

Figure Ml. 2 - THE DISCRETE-PLATE EQT ITIONS FOR A BINARY 
PLATE DISTILLATION COLUMN 



CHAPTER M2 



THE CONTINUOUS-SPATIAL EQUATION 

This chapter presents the steps in the procedure used to transform 
the discrete-plate equations of Figure Ml. 2 to the continuous-spatial 
equation (CSE) which is a continuous dynamic model of the concentration 
behavior of a binary plate distillation column. In general this con- 
version represents a significant step away from the quantitative accuracy 
of the discrete models because of the series expansions and approxima- 
tions involved. Usually, however, the solutions of the discrete models 
require comparatively large amounts of time to calculate, often on the 
order of half (about 15 minutes) of the major concentration time con- 
stant (B2.2 - Digital Computer). This makes such models of limited 
usefulness in control systems which utilize the model to predict 
transient response in order to correct for transients ahead of time. 

For these reasons a model is sought which requires less computation 
time to predict the dynamic response of the distillation column. 

The continuous-spatial equation (CSE) may be such a model. Ihe 
remainder of this thesis will be devoted to developing, investigating, 
approximating, and solving models of the continuous-spatial type. These 
models are not expected to behave quantitatively as near to the actual 
column response as do the discrete models, but it is anticipated that 
some sacrifice in numerical accuracy will result in a great savings in 
computation time and therefore better prediction for use in control. 

Ihe literature pertaining to continuous models is somewhat limited. 
The best single presentation on this subject in the literature of 
Chapter B1 is by Rosenbrock (R-12). It would seem that the most general 
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form of continuous model solved analytically in the literature is the 
Taylor Diffusion Model (Chapter IA) which results from approximations 
which, for the case of distillation, have been shown not to represent 
accurately the behavior of the column (G-3) (R-7). The most general 
form of continuous model solved numerically in the literature is the 
polynomial -class nonlinear model (Chapter M3). This model is solved 
by numerical methods which treat the continuous variable as a set of 
discrete points, and therefore the solution takes nearly as much time 
as solving the discrete equations (B2.2 - Transient or Time Analysis). 
It can be expected that most of the time savings resulting from solving 
continuous models will be a direct result of extensive analytical 
analysis and judicious approximation. This will be the objective of 
the following chapters in this thesis, 

The development of the continuous-spatial equation begins with 
the representation of the plate number n by a continuous-spatial 
variable x given in equation M2.1. The variable x is continuous from 

x = gn M2.1 

Where: n ** plate number 

g = plate spacing (ft. ) 

bottom to top of the column. The general discrete-plate equation is 
represented here from Chapter Ml as equation M2. 2, The transformation 

R d-U n = Vf (u n _i/ - Vf v.u n ) + Lu n+1 - Lu n M2, 2 
dt 

Where: u n = u R (t) 

from plate number to spatial variable means that u^ is redefined as 



- 88 - 



in equation M2. 3 and M2. 4, Using these functions, equation M2. 2 now 
becomes the partial differential equation of M2, 5. 



U n^ " u ( x,t ) M2.3 

u n +i( t ) ° u(x+g,t) M2. 4 

H 3u^x,t) = Vf [u(x-g,t)] - Vf [u(x,t)] M2. 5 

+ Lu(x+g,t) - Lu(x,t) 



The continuous-spatial equation M2. 5 is as exact as the discrete- 
plate equation M2. 2, only it is in a different form. Equation M2. 5 is 
too general for any analytical solution, thus, the terms in x-g and 
xtg are each expanded in a Taylor Series (G-3) , and up to the second 
order terms are retained as in equations M2, 6 and M2. 7. If these are 
now substituted into equation M2. 5 the result is given by equation M2, 8, 



Vf [u(x-g,t)] * Vf [u(x,t)] - gV d_ f [u(x,t)] 

2 X 

+ g 2 v d 3 

2 jjx* f [u(x,t)] M2. 6 



Lu(x+g,t) - Lu(x,t) + gL 3f :u(x,t)] 

3 x 

+ gfL^ 3 u(x,t) M2, 7 

2 fr? 



H 9u 

at 




[Lu + Vf(u)] + g [Lu 

3 X 



Vf(u)] 
M2. 8 



Equation M2. 8 can be written in a somewhat neater form with the 
use of several variable changes. Dividing through M2, 8 by V 
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and considering t in all of the above equations to be t = t., where 

c 

t is the column time variable, then the transformation to a new time 
c 

variable t is given by M2.9« Similarly the continuous spatial variable 

t = V t M2. 9 

H ° 

of the column is considered to be x c , and the change of variable on 
x is given by equation M2, 10, If the condenser and the reboiler are 

x = x /g M2, 10 

c 

now considered as the O'th and N + l'st plates respectively, then the 
ranges of the variables can be scaled as in M2. 11 and M2. 12, with the 
column height represented by (ft.) in M2. 13 for equally spaced plates. 



0 < x < C, 

— c — h 


M2. 11 


0 < x < 1.0 


M2. 12 


C h = (N+l)g 


M2. 13 



If all of these operations are applied to equation M2. 8, the result 
is equation M2.14. 

3u = 1_ 3 s |j3u + f( u ) j + ~ O “ f(u)] M2.14 

3t 2 3x‘ “ %'A 

Ihe complete model of the distillation column composition dynamic 
behavior using continuous spatial equations (CSE) can now be presented 
as Figure M2.1. This chapter has presented the steps in the procedure 
used to convert the DPE model of Chapter Ml to the CSE model of Figure 
M2.1. The CSE model is a highly nonlinear partial differential equa- 
tion of second order with nonlinear boundary conditions and as such 
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requires that approximations be made before any analytical results 
can be expected. Linear approximations to this model are the main 
subjects of the remainder of this thesis. The next chapter describes 
possible solution techniques applicable to the CSE and develops the 
Linear Polynomial -Coefficient Model (LPCM). 



Model Equation 




3u = 1 3 8 QBu + f(u)l + B ^Bu - f(u)] 
£t 2 $x 2 $x 




Where s = L u / V Upper x f < x < 1.0 




B^ = L^/V Lower 0 < x < x^. 




x^ = feed tray location 

Equilibrium 




f M - .gu . 

l+(a-i)u 




Boundary Conditions 




flu = f (u) - u x = 1.0 

fl* 


Top 


B flu - flf(u) = F (u-u f ) x = x f 

}x V 


Feed 


For: B = L /V and B., = L. /V 

u u' 1 1' 




B flu = f(u) - u x - 0,0 

$x 


Bottom 


Figure M2.1 - THE CONTINUOUS-SPATIAL EQUATION (CSE) 


DYNAMIC MODEL 
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CHAPTER M3 



SOLUTION TECHNIQUES FOR THE CONTINUOUS-SPATIAL EQUATION (CSE) 

This chapter presents a discussion of possible solution techniques 
applicable to the CSE, with particular emphasis on a nonlinear approx- 
imation to the LSE. The basic equation of the CSE is here repeated 
for convenience as equation M3.1. The CSE is a second-order, nonlinear, 
partial -differential boundary-value problem with nonlinear boundary 
conditions. An outline of possible techniques for solving the CSE 
is presented in Figure M3.1. In this chapter Parts A, B, C, and E 
of Figure M3.1 are discussed very briefly, and Part D is examined with 
the aim of developing the Linear Polynomial -Coefficient Model (LPCM) , 
which is examined in detail in Section 3(l)» 

3u = 1 3 2 [Bu + f(u)] + [Bu - f(u)] M3.1 

3t 2 fix 5 0X 

As stated previously, any distillation column model would consist 
of two separate CSE's of the form M3«l: one for the rectification or 

upper section of the column and one for the stripping or lower section. 
The two separate CSE 1 s are connected by the boundary condition at the 
feed tray, which represents the dividing point between the lower and 
upper sections of the column. Any complete solution for the dynamic 
composition behavior of a distillation column would therefore involve 
combining the solutions from two CSE's and their appropriate boundary 
conditions* The central problem to be considered now is that of 
solving the CSE* 
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In order for the reader to get some idea of how complicated the 
CSE really is, the equation is presented in expanded form as equation 
M3. 2, Considering Part E of Figure M3.1, the general theory of 
partial differential equations (B2.4 - Partial Differential Equation 
Theory) provides little help in finding the solution to the CSE. 

The general theory presents extensive and valuable results pertaining 
to existence of solutions, characteristic analysis, and methods of 
solution for linear, semi -linear, and quasi-linear equations, none of 
which can be applied to the completely nonlinear CSE. 



9u «* 1 [B + 0f (u)~ l 0 2 u + Ql ~o l2 f (u) , 0u 

5t 2 3 U 2 3x M3. 2 



- Bf (u) + b"| 9u 

0u 0x 



Where: u = u(x,t) 



B = constant = L/V 



f(u) 



3f(u) 

0 2 f (u) 



-US r 

+ (a-l) 



U 






[1 + (a-l)uj 
[1 + (a-l)uj 



2 

3 



Concerning Part A of Figure M3.1, it would seem to be equally 
difficult to simulate the entire set of discrete-plate equations as to 
simulate the CSE by any of the three methods listed. In fact, since 
the CSE is a partial differential equation, use of digital, hybrid, 
or analog equipment to simulate it requires that the CSE be broken up 
into discrete-spatial parts. Some time savings might result in appli- 
cation of an analog computer for the nonlinear spatial portion evaluated 



at discrete times. 
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A. Computer Simulation of the CSE 

1. Digital Simulation 

2. Hybrid Simulation 

3. Analog Simulation 

B. Application of Transformations to the CSE 

1. General Transform Theory 

a. Frequency Transforms; Laplace, Fourier, etc, 

b. Integral Transforms 

2, Change of Variables in the CSE to produce: 

a. More Easily Solved Nonlinear Equation 

b. Linear Equation 

c. More Easily Approximated Equation 

C. Application of Some Results in Pure Mathematics to the CSE 

1, Contraction Mappings and The Fixed Point Theorem 

2. Techniques for Finding the Fixed Point 

D. Approximation Techniques for the CSE 

1, Nonlinear Polynomial Approximations 

2. Linearization and Linear Polynomial -Coefficient 
Models (L?CM) of the CSE 

E. Application of the Theory of Partial Differential 
Equations to the CSE 

Figure M3.1 - POSSIBLE SOLUTION TECHNIQUES FOR THE CSE MODEL 



Solution of the CSE by application of some form of transforma- 
tion to the variables, Part B in Figure M3»l, is one area which might 
offer rewards with further investigation. Of course, frequency 
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transforms such as Laplace or Fourier can be eliminated from consider- 
ation immediately in the case of the general CSE because they depend 
upon linearity properties. However, these techniques might be effec- 
tively applied to the linearized CSE, and in fact, they are used for 
analytical and numerical solutions in the literature (B2.2 - Continuous 
Spatial Equations - Frequency Analysis.) 

In the application of a general transformation, an equation of the 
form M3. 3 is sought which will make the CSE a more easily solved non- 
linear equation, a linear equation, or a more easily approximated 
nonlinear equation. In addition, a variable transformation of the form 

F(s) = / f(u)K(s,u)du M3. 3 

o 

of equation M3. 4 could be applied to the CSE. The resulting complete 
set of possible transformations on dependent variables is then given 
by M3. 5. Transformations of this type are often used successfully for 

u = h(v) M3. 4 

F(s) = / ffh(v)"j K(s,v)dv M3. 5 

o 

solving problems in fluid mechanics and heat transfer. A search of 
these areas for transformations applicable to the CSE or some nonlinear 
approximation to the CSE might prove rewarding. 

Part C of Figure M3.1 presents an area of analysis which may 
possibly apply to a theoretical examination of the CSE. The major 
consideration here is the use of the Fixed Point Theorem of Brower 
(S-2l) (H-ll) to prove such properties as existence, uniqueness, and 
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continuity of solutions to the CSE. A description of the Fixed Point 
Theorem depends upon the concept of a contraction mapping. The appli- 
cation of a contraction mapping to the CSE depends upon considering 
the CSE as a nonlinear differential equation of the form of equation 
M3. 6 at each instant of time. 

du = f(x,u) M3.6 

dx 

The next step in applying the Fixed Point Theorem is the conver- 
sion of M3. 6 to an integral equation (B2.4 - Integral Equation Theory) 
presented in equation M3. 7. Then the right hand side of M3. 7 can be 

u(x) = u(x o ) + f[z,u(z) jdz M3. 7 

o 

considered as a general functional transformation (mapping of 
functions) A(v) given by equation M3. 8. Equation M3. 7 can next be 
written as a functional mapping of the function v onto the function 
u as in M3. 9. 



A(v) - v q + f(z, v)dz 

x o 


M3. 8 


A(v) = u 


M3. 9 



Now, any functional transformation A(u) which satisfies equation 
M3. 10 for any two functions u^(x) and u 2 (x) is defined to be a contrac- 
tion mapping. The Fixed Point Theorem of Brower (1912) then guarantees 

1*0^) - A(u 2 ) I < M I u q “ u 2 I M3. 10 

Where: 0 < M < 1 

for any contraction mapping that there exists a function u, the fixed 



-9 6 - 



point in the mapping, such that M3. 11 is satisfied. In other words 
the function transforms to itself. Equation M3. 11 is equivalent to 

A(u) = u M3. 11 

equation M3. 7, and thus the Fixed Point Theorem proves existence of 
a solution to this special case of the CSE if the CSE when written as 
a transformation can be shown to be a contraction mapping. 

The answer to the question of whether the CSE represents a con- 
traction mapping or not could only be determined by further detailed 
investigations of the CSE, However, since in the original column the 
variable u represented the concentration of the lighter component in 
the liquid, it could certainly be expected that u(x,t) would always 
remain in the range of M3. 12. Any distillation column model which 
has negative or greater-than-unity concentration solutions cannot be 
valid. In addition, all of the analytical solutions to approximate 
versions of the CSE in Chapter L4 satisfy equation M3. 12. 

0 < u(x,t) < 1.0 M3. 12 

The Fixed Point Theorem offers valuable information concerning 
the existence of a solution to a general nonlinear differential equa- 
tion, but it says nothing about the techniques for finding that solu- 
tion. It is expected that general solution techniques for an equation 
involving transformations and contraction mappings would have to be 
found in the mathematics literature (general topology) and applied 
to the CSE as a special case. It is highly unlikely that any such 
techniques could be applied to the CSE without some simplifying 
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assumptions or approximations being made. 

Thus, the analysis of the CSE now proceeds to Part D of Figure 
M3.1, beginning with nonlinear polynomial approximations to the CSE. 

The essence of nonlinear approximations lies in the assumption that 
the equations generated by approximating the coefficients of equation 
M3. 2 using the initial steady-state distribution will be valid for the 
transient behavior of the column. The validity of this assumption 
can only be determined by solving specific examples by numerical 
methods and comparing them to the CSE numerical solutions. 

One possible nonlinear approximation technique would be the use of 
n-th degree polynomials to represent the coefficients of equation M3. 2, 
The steady-state equation would then be a polynomial class (D-5, Ch.8) 
nonlinear ordinary differential equation of the form of equation M3. 13. 
The resulting nonlinear approximate model of the CSE would be given 
by equation M3. 14, where the and q^ are determined by approximating 
the coefficient functions of u and the equilibrium relationship by 
polynomials of sufficient degree for desired accuracy. 

A(u) dfu + B(u) du + C(u) / dm 2 + D(u) = 0 M3.13 
dx 2 dx ' dx' 

f * m i 

Where: A(u) = 2 A, u 

i=o 1 

B(u) = 2 Bn 1 

1=0 1 

C(u) = £ C.u 1 

i=o x 

y \ q _ i 

D(u) = 2 Du 

i=o 1 

A^, ®i* *^i» 3X6 functions of x 
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0 



M3. 14 



au - [p^u)] + 1 _ [p (u)] 

3t P ^ 

. . n i 

Where: P, (u) = £ p.u 1 

1 i=o 1 

, x m A 

P 2^ " V 

The nonlinear approximated version of the CSE presented in M3. 14 
is an analytically unsol vable equation in general. That M3. 14 can 
be solved using less computation time than that necessary for the CSE 
is very doubtful. It may be that M3. 14 for specific cases might be 
more easily programmed on an analog computer, but further investiga- 
tion for specific cases would be required to determine so. 

This chapter has presented a brief discussion of a variety of 
different methods which could be applied to the CSE. The results of 
the investigation are as expected; in order to proceed with any 
analytical analysis of the CSE, the equation must be linearized. 

The development of linearization techniques is presented in the next 
chapter (M4). 
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CHAPTER M4 



LINEAR APPROXIMATIONS TO THE CONTINUOUS-SPATIAL EQUATION (CSE) 

This chapter presents the steps in the linearization procedure 
to convert the CSE to a linear partial differential equation with 
spatially varying coefficients. The spatially varying coefficients are 
then approximated by n-th order polynomials in x, forming the Linear 
Polynomial -Coefficient Model (LPCM). The most important characteristic 
of linearized models is that the variable u(x,t) in the linearized 
model represents the incremental distribution resulting from the 
linearization about steady-state operation. 

The CSE is written in equation M4.1 in terms of u c (x,t), where 
u c (x,t) is the column concentration of the lighter component in the 
liquid. If M4.1 is linearized about the initial steady-state operation 
using equation M4.2, then the resulting linearized CSE model for u (x,t) 
representing the column deviation from the initial steady-state distribu- 
tion resulting from transient boundary condition variations is given by 
M4.3. The details of the linearization procedure are presented in 
Figure M4.1. Equation M4.3 can then be written as a linear partial 

Bu 1 B 3 fBu + f(u )"] + 3 [Bu -f(u c )l M4.1 

^=20? L C C ^ ° 

u c (x,t) = u^x) + u(x,t) M4.2 

differential equation of second-order with spatially varying coeffi- 
cients by expanding the terms as in equation M4.4. 

Bn = 1 3 s [jBu + m(x)u^ + [jBu - m(x)u^| M4.3 

0t 2 0*> 0x 



- 100 - 



Development of the Linearized CSE 



A (u) => 1 cLf_ [Bu + f(u)] + b_ [Bu - f(u)1 - 3u 

3 x $t 



A(u ) = A(u ) + 3 a 

u = u c - Ui 

A( u i ) = A(u c ) = 0 



(u) 



u. 

1 



dk 

a^ 

3a 



(u) » 0 



- 1 a 2 U (B + m(x))u^j + 3_ [(B - m(x))u[] - 3_u = 0 



2 3 x 



ax 



Development of the Boundary Conditions 
A^.(u) = ' d - 1) f(u) + u = 0 

g x 

Expanding similarly to above: 
ft - 1\ (m(x)u) + u = 0 

f dm(x) - m(x) + 1 u + m(x) du = 0 
dx gx 

A (u) = B [u - f(u)] - F (u - u f ) = 0 

3 X v 



at 



x = 1.0 Top 



x = 1.0 Top 



x = Xj. Feed 



Expanding as above: 

-[B dm(x) + Fj u + jjL - m(x)^ bu = -F u^(t) 
dx V a x v 





X = X« 


Feed 


Where: u^(t) = u^ c (t) - u^ 


Upper and Lower 


A b (u) = (B n t- + 1) u - f(u) 

a x 

Expanding as above: 


x = 0.0 


Bottom 


£l - m(x)”| u + 3, 3u = 0 


x = 0.0 


Bottom 



1 ~ 

a x 



Figure Mfr.l - DETAILS OF THE CSE LINEARIZATION 
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h.- i [b ♦.(*>! &. + [B - m(x) + dm(x)~| du 
2 3X 2 dx dx 

+ [jL d 3 m(x) - dm(x)~ | u M4.4 

2 dx 3 dx 

It is very Important to realize the nature of the function m(x). 
Ihe function m(x) is defined in equation M4.5 and represents the slope 



m(x) = ^f(u) M4.5 

0u 

u = u^x) 

of the equilibrium curve as a function of x. In equation M4.3 the 
constant B represents the slope of the operating lines on the McCabe- 
Thiele diagram of Figure 12.2 in the upper and lower regions. The 
function m(x) represents the slope of the equilibrium curve in those 
same regions. Almost all approximation techniques applicable to 
equation M4. 3 use these characteristics and assume that m(x) is given 
by an equation of the form of M4.6 or M4.8 (R-12) (G-3) (B2.2 - 
Analytical Analysis). With the assumption of equation M4.6 equation 
M4.3 becomes equation M4.7 which is the Heat Equation. 



m(x) = constant = B M4.6 

9u = B 9fu M4. 7 

0t 0x 2 

It is rather obvious that assumption M4.6 reduces the linearized 
equation M4. 3 to the simplest diffusion equation, thereby neglecting 
most of the behavior which is unique tc distillation. Distillation 
is a diffusion process, and it is somewhat reassuring to find that 
the model M4. 3 reduces to a diffusion model. However, it is certain 
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that a model of the form of K4.7 could not be used to accurately 
predict column behavior for the purpose of controlling the column. 
Equation M4.7 is solved analytically and numerically for B = 1.25 in 
Chapter L4 and represents the "bottom of the ladder" in any hierarchy 
of approximate models of the CSE. 

The next step upward in complexity results from assumption M4.8 
resulting in equation M4.3 becoming M4.9. Equation M4.9 is the Taylor 
Diffusion model equation and is discussed in detail in Chapter L4. 

This equation is still too approximate to use for any control applica- 
tions and so, a more sophisticated representation of equation M4.3 is 
sought. 

m(x) = constant / B M4.8 

Bn = P}3 3 u + Pg^u 

3 t ~W TT ^ 

Where: P. = (B + m)/2 
P 2 * (B - m) 

Suppose the spatially varying coefficient functions in equation 
M4.3 are approximated by n-th degree polynomials in x. Then the 
resulting complete model would be the Linear Polynomial -Coefficient 
Model (LPCM) presented in this chapter as Figure M4.2 and in 
Chapter LI as Figure Ll.l. The boundary function constants in 
Figure Ll.l are equivalent to the coefficient functions, evaluated at 
the boundaries, in Figure M4.2. 

The central purpose of this thesis is to suggest that the LPCM of 

Figure M4.2 can be used as a distillation column dynamic model of 

sufficient accuracy and rapid solution as to be useful in controlling 
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the column and to present an analytical solution technique for it. 
This is the subject of Section 3(b). 



Model Equations 

£u = 3 s pP-,(x)u'] + [P ? (x)ul 

at a * 3 a x 

Where: P.j,(x) = £ PjX 1 = [B + m(x)]/2 

p o( x ) = 2 q.x 1 = [B - m(x)*] 
c i=o 1 

Upper and Lower Equations, One for Each Section 
of the Column 

Boundary Conditions 

( jim(x) - m(x) + l] u + m(x) 9u = 0 x = 1.0 Top 
dx 0x 

-[B dm(x) + F] u + Tl - m(x) j “ ^E u„(t) 
dx V 3x V * 

x = Feed 

(Upper and Lower) 

pL - m(x) j u + B.0u =0 x = 0,0 Bottom 

0x 



Fifuire M4.2 - THE COMPLETE LINEAR P0LYN0MIAL- 
COEFFICIENT MODEL (LPCmT 
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SECTION 3 



THE LINEAR POLYNOMIAL - COEFFICIENT MODEL (L) 

LI ANALYTICAL SOLUTION OF THE LPCM BY INTEGRAL EQUATION TECHNIQUES 
L2 THE PARTIAL LINEAR POLYNOMIAL BOUNDARY VALUE PROGRAM (PLPBV) 

L3 DETERMINATION OF THE LPCM FOR EXAMPLE DISTILLATION COLUMNS 
L4 ANALYTICAL SOLUTIONS TO APPROXIMATED COLUMN EQUATIONS 
L5 OTHER SOLUTION TECHNIQUES APPLICABLE TO THE LPCM 



"DAZZLED BY THEIR ABILITY TO DO ELEMENTARY THINGS AT TREMENDOUS 
SPEEDS AND TO PUT THESE TOGETHER IN STRUCTURES OF DAUNTING COM- 
PLEXITY, SOME HAVE ALLOWED THE TERM 'GIANT BRAINS' TO GAIN 
CURRENCY AND, SEDUCED BY THE SIREN SONG OF SO SENSELESS A 
SOBRIQUET, HAVE SURRENDERED THEIR BIRTHRIGHT OF RATIONAL THOUGHT 
FOR A POTTAGE OF PUNCHED CARDS." - R. ARIS (A-13) 



THE CENTRAL PURPOSE OF THIS THESIS IS THE DEVELOPMENT, 
PRESENTATION, SUGGESTED SOLUTION TECHNIQUE, AND EVALUATION OF 
THE LINEAR POLYNOMIAL - COEFFICIENT MODEL (LPCM) OF THE DYNAMIC 
BEHAVIOR OF A BINARY PLATE DISTILLATION COLUMN. THE FUNCTION 
OF THIS SECTION IS THE PRESENTATION AND UTILIZATION OF THIS MODEL 
AND THE DEVELOPMENT OF AN INTEGRAL EQUATION SOLUTION TECHNIQUE 
FOR IT. 
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CHAPTER LI 



ANALYTICAL SOLUTION OF THE LPCM BY INTEGRAL EQUATION TECHNIQUES 

Ihe purpose of this chapter is to present a somewhat involved 
analytical solution of the Linear Polynomial -Coefficient Model (LPCM) 
defined by equation Ll.l with boundary conditions LI. 2, LI. 3* and 
LI. 4. The LPCM is a linear, second-order, parabolic, partial- 
differential, two-point, boundary-value problem and, as such, cannot 
be solved analytically or numerically with complete generality (B-8). 
The best that car be done analytically is to develop a solution 
technique which allows reasonable assumptions and approximations to 
greatly simplify the analytical manipulations leading to a simplified 
solution. The best that can be done numerically is to develop and 
employ programs specifically related to specialized cases of the LPCM 
and then to utilize those programs on digital, hybrid, or analog 
equipment to compute a simplified solution. In either case, the 
solution may or may not truely represent the behavior of the distilla- 
tion column from which the model was developed. 

Equations of the form of the LPCM are usually solved numerically 
by a variety of different techniques. Numerical solutions to the 
LPCM are described in Chapter L5 and a survey of the literature in 
the Bibliography, Chapter Bl, pertinent to solving the LPCM is 
presented there. No attempts have been made to solve the LPCM in a 
purely numerical fashion in this thesis, although the suggestion can 
be made that such attempts might prove to be exceedingly profitable 
in terms of greatly reduced solution computation time when using the 

LPCM as part of a column control system. 
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The solution technique to be described in this chapter is a 
sequence of analytical manipulations which lead to am expression 
which is then to be evaluated numerically. The technique begins with 
expressing the response in terms of steady state and transient 
portions. The technique of separation of variables is then employed 
for the tramsient partial differential equation. The spatial 
differential equation portion resulting from the separation of 
variables is then solved by converting it to a Liouville Normal- 
Form equation and then to a homogeneous Fredholm II-integral equation. 
The total solution is then expressed in terms of the eigenvalues and 
eigenfunctions resulting from the solution of the integral equation. 

The anticipated advantage occuring from these analytical manipulations 
is that the final numerical evaluation in this technique may take 
much less computation time than solving the LPCM purely numerically. 

Ll.l TRANSIENT AND STEADY-STATE EQUATIONS FROM THE LPCM 

The Linear Polynomial -Coefficient Model is defined by equation Ll.l 
with boundary conditions LI. 2 , L1.3» and LI. 4. This model represents 
the deviation from steady state of the concentration of the lighter 
component in the distillation column subject to step changes in the 
end-point compositions. The complete presentation of the development 
of this model has been given in Section 2(M). 

The first step to be taken in the analytical solution of Ll.l 
is the realization that, at any given instant of time, there are 
three separate spatial concentration distributions implied by the 
model. First of all, there is the initial steady state distribution 
u^(x) in the column prior to the step inputs. Since the LPCM 
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Ll.l 



IfL Di(x)u] + 3 _£p 2 (x)u] - 8u 

7 ** $x Qt 

Where* u - u(x,t) 

P x (x) - 8 p x 3 - 
A i»o 1 

p 2 (x) 

A 11 u(b,t) + A 12 u x (b,t) . A^ + A^U^t) LI. 2 

A 2 i u ( a i"t) + A 22 u x (a,t) «# A^3 + A 2 ^U_^(t) LI, 3 

A 11 A 22 ” A 12 A 21 ^ ® LI, 4 

Figure Li.l - THE LPCM 



represents an approximation to an equation which has been linearized 
about an initial steady state, this portion should be zero if the 
model is to be valid, but cannot be assumed zero for any general 
model. Thus, u^(x) must satisfy equation LI. 5 with boundary conditions 
given by LI, 6 and LI. 7, 

3 3 r P 1 (x)u 1 (x)1 + 3_ [P 2 (x)u i (x)] = 0 LI. 5 

9^ 3x 

+ A i2^i^ b ^ “ ^3 LI. 6 

A 2 l u i( a ) A 22 u i a A 23 LI. 7 

The second distribution implied by the model is the final steady 

state u s (x). This represents the lighter component concentration in 

the column after the step inputs and after all transients have had 
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time to die out. The final steady state u g (x) must also satisfy 
equation LI, 5* but with the boundary conditions given by LI. 8 and 
LI. 9. 

A^Cb) + A 12 i 3 (b) " A-^3 + L1.8 

A 21 u s(a) + A 22 u s^ a ^ “ A 23 + a 24 LI. 9 

The third concentration distribution is the transient response 
u(x,t) - (Note: the same notation is used for the transient 

response as was used in Ll.l for simplicity) - which must satisfy 
equation Ll.l, but with the boundary conditions LI. 10 and LI. 11 
and the initial condition LI. 12. In the interest of having a com- 
pletely defined problem and also from the practical constraint of 
proper model behavior, the final condition LI. 13 must also be 
included. 



A}ju(b,t) + A^u^b.t) 


=S 0 


LI. 10 


^21 u ' ,a *t) ^ AggU^ajt) 


=2 0 


LI. 11 


u(x,0) = u Q (x) - u i (x) 


- u a (x) 


LI. 12 


O 

II 

IT 

X 




LI. 13 



The net result of considering these three spatial distributions 
is that the LPCM solution now is seen to be the solution of three 
separate problems* two ordinary differential problems LI. 5 with 
different boundary conditions and one partial differential problem 
Ll.l with conditions LI. 10 - LI. 13. 
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The general usefulness of the LPCM must lie in the capability 
of using step response solutions to approximate the response to any 
arbitrary input functions of time. In utilizing the model in this 
general manner, the initial steady state u^(x) cannot be assumed to 
be zero and the model must be solved as the three separate problems 
described previously. 

For the purposes of presenting and solving the LPCM in terms of 
single step inputs, however, the assumption will be made in this 
thesis that and are both zero and, thus, that the initial 
steady state distribution u. (x) is zero throughout the column, i.e. 
that the initial deviation from steady state is zero. Solution of 
the general problem implies solving for u^(x) in LI. 5 and using it, 
along with u (x), in the transient boundary condition LI. 12. 

The assumption that u^(x) = 0 for the purpose of evaluating 
the model reduces the solution to one partial differential problem 
Ll.l in the transient response u(x,t) and one ordinary differential 
problem LI. 5 in the steady state response u (x) with boundary con- 
ditions given by LI , 14 and LI. 15* Once the transient and steady state 
solutions have been determined, the total model solution u (x,t) for 

HI 

A ll u s( b ) + A 12^ b ) " ^4 LI. 14 

A 2i u s ( a ) + A 22^ a ) “ A 24 LI. 15 

step inputs can be expressed as L1.16* Then u m (x,t) will satisfy 
the LPCM with * A^ - 0. The development of these solutions 
begins with the application of the technique of separation of 
variables to the transient partial differential equation Ll.l. 
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u m (x,t) = u s (x) + u(x,t) 

LI. 2 SEPARATION OF VARIABLES APPLIED TO THE LPCM 



LI. 16 



The method of separation of variables rests on the assumption 
that the solution u(x,t) can be separated into a product of two 
functions X(x) and T(t) as in eq. LI, 17. This assumption is the 

u(x,t) - X(x) • T(t) LI. 17 



single most universally used analytical technique for solving 
partial differential equations. If this expression for u(x,t) is 
then substituted into the LPCM, the resulting expression is 
equation LI. 18, where P^ and are defined in LI. 32 and LI. 33. 

? 1 XT + P^T + P^XT - X f £ LI. 18 

If this expression is divided by XT, then the resulting expression 
in equation LI. 19 represents a function of x equated to a function 
of t for all values of x and t. Thus, the two expressions must 
equal a constant - K 3 . 

The applicability of this technique to a given partial 
differential equation rest3 on the capability of separating the 



P^X + P 3 X + P^X _ T = _ k2 
X T 



LI. 19 



resulting equation, as in equation LI. 18, This step is, in general, 
not possible when the given partial differential equation is non- 
linear as are, for example, the continuous-spatial equation and 
polynoraial-class partial differential equations (Davis, D-5, p.213). 
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Ihe general condition for separability of a partial differential 
operator of the form of LI. 20 has been shown by Murray-Lasso (M-5) 
to be the commutation of the operators H and H, expressed in eq, 

LI. 21. In this case H x , 1^, and L^ t are partial differential 
operators in x, t, and both x and t, respectively. In the case 

L x,t “ K t [u(x,t)] + H x [u(x,t)] LI. 20 

H x C H tC u >l - H t [**(«)] LI. 21 



of the continuous-spatial equation H t and H x are defined by 
equations LI. 22 and LI. 23, where f(u,u , u ) is a nonlinear function. 

A XX 



H t 00 = 3h 




at 


L1.22 


H xM ■ f ( u 'V U xx) 


LI. 23 



That these two operators do not commute is shown in eq. LI. 24 and 
therefore the continuous-spatial equation is not separable. The 



3 [H v (u)l - 3 h v 3 

aT L x v '-I X 0U = u r, -i 

3* 3t * LU tJ 



LI. 24 



LPCM, however, does satisfy eq. LI. 21 and is separable, as has 
been shown in eq. LI. 19. 

The application of the separability assumption then results 
in the transformation of the partial differential equation in two . 
variables to a set of two ordinary differential equations, each 
of which depends upon the separation constant -K 3 . The two 
resulting ordinary differential equations are presented as the 
spatial equation LI. 25 and the time equation LI. 26. 
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LI. 25 



P-Jx + P 3 X + P^X - -K 3 X 

T - -K a T LI. 2 6 

LI. 3 SOLUTIONS OF THE TIME EQUATION 

The time equation LI. 26 is an eigenvalue problem whose eigen- 
function solutions are easily found. Since u(x,t) in equation Ll.l 
represents the transient portion of the desired solution, the 
separation constant has been chosen as -K 3 in order that each 
eigenfunction of the transient response approach zero for large 
values of time as in equation LI. 13. The eigenfunctions which satisfy 
equation LI. 26 are shown as functions of the eigenvalues -K^ 3 in 
equation LI. 27. 

T Q (t) - expC-K/t) LI. 27 

LI. 4 SOLUTIONS OP THE SPATIAL EQUATION 

The eigenvalue problem presented by the spatial equation LI. 25 
and its associated boundary conditions is much more difficult to 
solve than the time equation problem. In fact, statements about 
the existence, uniqueness, and continuity of solutions for this 
problem cannot be made in general (B-8). Statements of this nature 
require restrictions on the equation or the boundary conditions. 

There are some existence and uniqueness theorems for restricted 
cases presented in the literature (3-9) (K-l) (H-ll), however, the 
existence, continuity, and completeness of eigenfunction solutions 
to regular Sturm-Liouville problems will be assumed and utilized in 
this section (B-14) (l— 2) . 
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Several very simple cases of the spatial equation LI. 25 can 
be solved analytically. These simple cases are for P-^ P y and P^ 
set equal to constants, and examples of these are solved in Chapter L4. 
In all of these simple cases and, in fact, for the problem LI. 25 in 
general, the final steady state distribution in the transient solution 
boundary condition LI, 12 is expanded in terms of the eigenfunctions of 
problem LI. 25. For reasonably well behaved systems, such as a binary 
distillation column, it is expected that only the first few, say ten 
(ID) or less, of the eigenvalues are really crucial to the transient 
behavior of the distillation column. The reason that only the 
smallest eigenvalues -Kjj 3 are important to the total transient 
behavior is that the time-eigenfunctions in equation LI, 27 approach 
zero for very small values of time when the eigenvalues -K^ 3 become 
larger. As mentioned in Chapter L4, the analytical solution eval- 
uations of the simple cases of the LPCM always required less than 
ten (10) eigenfunctions for three-decimal -place accuracy. 

These considerations offer promise of greatly simplified 
solution expressions for the LPCM if an analytical technique can be 
developed which offers a solution in terms of the first few eigen- 
values and eigenfunctions. As mentioned in Chapter L5, there are 
numerous analytical and numerical procedures which offer this rep- 
resentation; however, conversion of LI. 25 to an integral equation and 
then solution of the integral equation either by separation of the 
kernel (Appendix A2) or by approximating the integral equation by a 
system of linear algebraic equations (Appendix A3) will be investigated 
in this thesis. 
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The integral equation technique consists of converting the 
differential equation to the Liouville Normal -Form equation and then 
converting the Liouville Normal -Form equation to a Fredholm II- 
integral equation. The first conversion procedure begins with a 
test for self-adjointness and a transformation of LI. 25 to a regular 
Sturm-Liouville equation LI. 40. The end result of these manipula- 
tions is the set of spatial equation eigenvalues K n in equation LI. 49 
and the set of eigenfunctions W n (z) in either equation A2.15 or 
A3. 10 with the transformation relations Ll,_44, LI. 45, LI. 37. and 
11,39 leading to the eigenfunction solutions X n (x). The following 
sections present the development of these transformations and the 
equation solution procedures, 

LI. 5 SELF-ADJOINT CONDITIONS AND TRANSFORMATIONS 

If the spatial equation LI. 25 can be shown to be self-adjoint, 
then the spatial eigenvalue problem is a regular Sturm-Liouville 
eigenvalue problem. If the spatial equation is not self-adjoint, 
then the variables can be transformed into a self-adjoint equation 
which then represents a regular Sturm-Liouville problem. In an 
analytical analysis such as employed in this chapter, it is' essential 
that at some step the equations be expressed in a format for which 
real eigenvalues and complete eigenfunctions are guaranteed. The 
regular Sturm-Liouville problem is such a format (B-14) (1-2). 

The necessary and sufficient condition that the spatial equation 
LI. 28 be self adjoint is expressed by equation LI. 29. When this 

P 1 (x/X + P 3 (x)X + [P 4 (x) + K«] X - 0 LI. 28 

Pl(x) - P^Cx) 
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L1.29 



condition is applied to the coefficient polynomials in equations 
LI. 30* LI. 31, LI. 32, and LI, 33* the polynomial coefficient relations 
LI. 34 and LI. 35 result. The self-adjoint form of equation LI. 28 
is then given by equation LI. 36 . 

P x (x) - 5*p, 1 x j ' 1 = £ p x 1 LI. 30 

x j=l J”- 1 - i=o 1 

, v n+1 -i,.! n i 

p 2 (x) - 5-i^j-i x - g, 0 v u -3i 

P 3 (x) = q 0 + ^ =1 (2j Pj + q^x* 5 ” 1 LI. 32 

Pq(x) * q x + |^(j+l) [j Pj+1 + q^x^’ 1 LI. 33 

q j-l + JPj “ ° J " 1* n LI. 34 

q n - 0 LI. 35 

d_ rp-,(x)x] + P k (x)x - -x 2 x Li .36 

dx L 4 

If the spatial equation LI. 28 turns out to be non-self-ad joint, 
then the transformation relations LI. 37* LI. 36* and LI. 39 (M-2) 
result in the self-adjoint equation LI, 40, This transformed 



P(x) - exp P^Cx) dx] 

a piSt 


LI. 37 


Q(x) - Pq,(x) # P(x) 

h &) 


LI. 38 


R(x) « P(x) 


LI. 39 


d_ [P(x)X] + Q(x)X = -K 3 R(x)X 


L1.40 



dx 



equation LI. 40 is then seen to be of the same form (Sturra- 

Liouville) as the self-adjoint equation LI, 36 . The spatial 

-II6- 



eigenvalue problem thus reduces to a regular Sturm-Llouville problem 
with the differential equation expressed by LI. 40 and the boundary 
conditions by L1.41, LI, 42, and LI. 43. 



A-j^XCb) + A 12 X(b) - 0 
A^xCa) + A 22 X(a) - 0 

A 11 A 22 “ A 12 A 22 ^ ° 



LI. 43 



LI. 42 



LI. 41 



The next step in the solution procedure is to convert the 
regular Sturm-Liouville equation LI. 40 to the Liouville Normal-Form 
equation LI. 46. Equation LI, 25 could have been converted directly 
to an integral equation but this was not done for several reasons. 
Firstly, the analytical manipulations in the conversion procedure 
did not lead themselves to the approximation techniques so necessary 
for finding the final solution. Secondly, and much more importantly, 
there is no guarantee that a direct transformation would produce an 
integral equation with a continuous kernel. In fact, in most cases 
it would not. If either the kernel separation method in Appendix A2 
or the linear algebraic method of Appendix A3 were employed for 
solving an integral equation with a discontinuous kernel, there would 
be no guarantee that the eigenvalues would be real and distinct. Such 
a guarantee can be given if the kernel is continuous and symmetric. 

LI. 6 TRANSFORMATION TO THE LIOUVILLE NORMAL-FORM EQUATION 

Transformation on both the dependent and independent variables 
of equation LI. 40 can be used to greatly simplify the regular Sturm- 
Liouville equation (B-14) (l-2). If new variables W and z are defined 
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as in equations LI, 44 and LI. 45, then the equation LI. 40 is transformed 
to that of equations LI. 46 and LI. 47, where the primes are derivatives 



V(*)R(xT 



LI. 44 
LI. 45 



with respect to x and the equations are expressed as functions of z. 



d 3 W 

dz* + [X - M(z)]W = 0 LI. 46 

a i B n , 12 

«(z) - P. [(£) ♦ («-) ♦ 1 (f) 

+ 1 (£■)$■) - i(2i) 2 ] + $. Li* i*7 

2 F K 4 K R 

Whereas the original equation LI. 40 was defined over the interval 
a < x < b, the transformed equation LI. 46 is valid over the interval 
0 < z < c, where c is defined in equation LI. 48 and X in equation LI. 49. 

c “ J'a Vtjgp* LI. 48 

X - K 2 LI. 49 



The boundary conditions of LI. 41 and Li. 42 are also transformed to 
the new boundary conditions on W(z) expressed in equation LI. 50 



through LI. 56, where the dot is d/dz and the prime is d/dx 


D n W(c) + D 12 W(c) - 0 


LI. 50 


d 21 w(o) + d 22 w(o) - 0 


LI. 51 




LI. 52 
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LI. 53 



°12 " A 12 




D 21 " A 21 ~ 




LI. 54 




LI. 55 



D 11 D 22 " D 12 I) 21 ^ 0 



LI. 56 



LI. 7 SOLUTION OF THE LIOUVILLE NORMAL-FORM PROBLEM 

Now that the spatial differential equation LI. 28 has been 
transformed to the Liouville Normal-Form, the next step in the 
solution of the LPCM consists of finding the eigenvalues and the 
eigenfunctions W i (z) for equation LI. 46 with the boundary conditions 
LI. 50 stnd LI. 51. As discussed in subsection LI, 4 the technique to 
be used is that of converting LI. 46 to an integral equation and then 
solving the integral equation. 

The conversion procedure is presented in Appendix A1 and two 
alternative solution procedures are presented in Appendices A2 and 
A3. The separable kernel procedure presented in Appendix A2 could 
be employed for specialized cases of the LPCM. The solution procedure 
in that case would be to approximate the kernel A1.17 by a finite 
series A2.2 and then to apply the techniques of Appendix A2, This 
would involve a great deal more analytical analysis for any given 
LPCM, but it is anticipated that the resulting solutions would be 
more accurate than those obtained by the methods of Appendix A3. 

The relative merits of these two solution procedures in terms of 
solution accuracy and computation time would have to be examined for 
any given LPCM. 
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The second integral equation solution procedure in Appendix A3 
consists, in essence, of approximating the integral by a finite sum. 

Ihe procedure is easily programmed on a digital computer and is a 
very direct method. For these reasons, the procedure in Appendix A3 
will be used to evaluate the LPCM in this thesis. 

Once the eigenvalues and eigenfunctions X^x) of the spatial 
equation have been found, the next step in finding the total solution 
to the LPCM is to ensure that the initial condition LI. 12 is met. 

This requires that the initial distribution u Q (x) be expressed in 
terms of those eigenfunctions X n (x). 

LI. 8 THE EIGENFUNCTION EXPANSION 

The eigenfunction expansion begins by expressing u Q (x) in equation 
LI, 12 as a sum of constants c^ multiplied by each eigenfunction X^(x) 
as in equation LI. 57, when n is the number of eigenfunctions. Since 

u o( x ) " | = - c i X i(x) L1 *-57 

the integral equation solution technique only gives n eigenfunctions, 

then equation LI. 57 can only be satisfied for n points x. in the 

J 

interval a < x < b. Thus, equation LI. 57 must be discretized to 
equation LI. 58. Equation LI. 58 then represents n-simultaneous equations 

u o(*j> ■ L1 - 58 

in n-unknowns which can be expressed as the matrix equation LI. 59. 

The solution to LI. 59 is then given by L1.60. Once the matrix C has 

U Q - XC LI. 59 

C - x“ 1 u 

o 
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L1.60 



been found, the total solution can be expressed in terms of the 
transient and steady state portions, 

LI. 9 THE TOTAL SOLUTION TO THE LPCM 

The total solution to the LPCM can now be expressed in terms of 
known functions u s (x) and u(x,t). The steady state problem LI. 5 
with boundary conditions LI. 14 and LI. 15 is a standard linear boundary 
value problem for which solutions and, in fact, computer subroutines 
exist. The steady state solution is determined in this thesis by 
Scientific Subroutine Program - LBVP in reference (l-3) and is 
described in greater detail in Chapter L2. 

The major portion of the LPCM, for which computer programs have 
not been developed, was the solution to the transient boundary value 
problem. For this reason, this portion of the LPCM received the 
greatest portion of solution effort in this thesis. It is expected 
that the steady state equation could be solved by the same integral 
equation techniques (leading to a final matrix inversion rather than 
an eigenvalue problem) as were used for the spatial equation portion 
of the transient response. The conversion would be more involved 
because the boundary condition constants and A ^ in equations 
LI. 14 and LI, 15 are not zero. The conversion of both the steady state 
equation and the spatial equation of the transient response, and the 
incorporation of the solution technique for both of these problems 
in the same program would undoubtedly result in a great savings of 
computation time over that spent using the LBVP subroutine and the 
integral, equation techniques separately. This could only be determined 
from further studies. 
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The efforts of the previous sections of this chapter have led 
to the transient solution to the LPCM due to step inputs. The total 
solution was given in equation Ll.l6j the transient portion of that 
solution can now be given as equation L1.61. The total solution to 

u(x,t) - £ c.X (x) exp(-K. s t) L1.61 

i=i 1 ± 

the LPCM is then given by equation LI. 62, 

u^x.t) * u s (x) + | CjX^x) expC-KjH) LI. 62 

The next chapter (L2) develops a digital computer program for 
performing these manipulations and calculating the solution to a 
special version of the LPCM. 
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CHAPTER L2 



THE PARTIAL LINEAR-POLYNOMIAL BOUNDARY VALUE PROGRAM (PLPBV) 

The very general solution technique developed for the LPCM in 
Chapter LI is based upon mathematical transformations whose analytical 
expressions, for general n th order polynomial coefficients, are too 
complicated to use in a numerical solution evaluation. This chapter 
presents an application of the integral equation technique of 
Chapter LI to the first-degree polynomial model (PLPBV) and a digital 
computer program for evaluation of the solution. 

L2.1 PRESENTATION OF THE MODEL ( PLPBV ) 

The Partial Linear-Polynomial Boundary Value (PLPBV) model is 
presented in Figure L4.X. This model is a subcase of the LPCM of 
Figure Ll.l in which P£(x) is a linear polynomial and P^(x) is a 
constant. This model represents one step upward from the simplified 
LPCM examples of Chapter L4 and, at the same time, the simplest form 
of the LPCM which still has a spatial varying coefficient. 



31 Pi Mu] 


' * ^ L P 2^ x ^ u J “ — 


L2.1 




ax 




Where: 


u * u(x,t) 






P x (x) « P 0 ; p Q > o 






p 2(x) - q Q + qi* 




AjMi(b,t) + 


Ai2U x (b,t) = A^^U^(t) 


L2.2 


A 21 u(a,t) + 


A 2 2 u x ( a » t ) “ A 24 U -1^^ 


L2.3 


Figure L2.1 - THE PLPBV MODEL 
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The development of the solution to this model now parallels that of 
Chapter LI for the general LPCM except that the transformation 
relations can he expressed specifically in terms of the polynomial 
coefficients. 

L2.2 THE TRANSIENT AND STEADY-STATE EQUATION FROM PLPBV 

The first step in solving the PLPBV model is the expression of 
the separate transient and steady-state problems as in subsection Ll.l 
sind equation L1.16. The transient and steady-state problems are 
presented in detail in Figure L2.2, where u^Xjt) is the total solution 
to PLPBV. 



u ( X , t ) - u (x) + u(x,t) 
m s 


L2.4 


Transient Portions 




P.u + P~u + P,,u « u. 


L2.5 


1 xx 3 x 4 t 


where: P-, =° p 

x 0 




P 3 = q 0 + 
P 4” *1 




A^uCb.t) + A i2 u x( b,t > " 0 
A 21 u(a,t) + A 22 u x (a,t) » 0 
u(x,o) * -u g (x) 
u(x.oo) = 0 




Steady-State Portions 




p i“s + P 3“s * Vs " 0 

where s P, = P 
1 0 


L2.6 


P 3 = q Q + qjX 

p 4 = q l / 

(Continued 


on next page) 
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^ll u s^ + A 12 u s^ b ) " A 14 
A 2 l u s (a) + A 22 u s (a) - A 24 

Figure L2.2 - TRANSIENT AND STEADY-STATE 
PROBLEMS OF PLPBV 



L2.3 SOLUTION OF THE STEADY STATE PROBLEM 

The solution of the general LPCM steady state has been discussed 
in subsection LI. 9. The subroutine LBVP, Linear Boundary-Value 
Problem, and the associated subroutines applicable directly to the 
general LPCM are presented in Appendix A5 for completeness and for 
the convenience of the reader. Subroutines LBVP and GELG are taken, 
less comment cards, directly from the literature (l-3)» whereas 
subroutines AFCT, DFCT, and FCT have been written specifically for 
the general LPCM. The user must supply his own output subroutine 
OUTP as per reference (l-3). 

Once LBVP is called in the main portion of PLPBV then it calculates 

the steady-state values u (x.) for use in equations LI, 58, for the 

s J 

final transient eigenfunction expansion, and In LI. 62, for the total 

solution evaluation at the discrete points x.. The user OUTP sub- 

routine must be set up to return u (x.) values to PLPBV at exactly 

s' y 

the desired x^. 

Some comments concerning the efficiency of LBVP in the applica- 
tion must be made. It is obvious that since both FCT and DFCT return 
zero values that a simpler version of LBVP for this special case 
could be written a r possibly found in the existing literature by 
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further searching. As mentioned in LI. 9, the integral equation 
techniques might also be applied here. 

For the initial purposes of evaluating the LPCM and testing the 
transient portions of PLPBV, the steady state constants were fed in 
as Column 9» Matrix 3, in Appendix A4, Subroutines LBVP, AFCT, DFCT, 
and FCT have all been compiled, tested, and shown to operate properly. 
L2.4 SOLUTION OF THE TRANSIENT PROBLEM; ANALYTICAL ANALYSIS 
Ihe solution to the transient portion of PLPBV defined in 
Figure L2.2 will be obtained by applying the transformation rela- 
tionships to L2.5 to convert the spatial equation equivalent to LI. 25 
to the Fredholm II-Integrai Equation A1.21, This analysis will be 
divided into an analytical analysis in which the specifics of the 
transformations will be presented and a computational analysis in 
which subroutine PLPBV will be explained and presented. 

The analytical analysis of the transient portion of PLPBV begins 
with separation of variables. The time equation eigenfunctions 
resulting from the application of the separation-of -variables 
technique are the same in PLPBV as they were for the general LPCM 
in equation LI. 27, The resulting spatial equation is given below as 
equation L2.7. 

P,X + P^X + [Pju, + K 2 ] X = 0 L2.7 

where: P. * p Q 

? 3 = i 0 + 

P 4 = *1 
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The next step in the analysis of the spatial differential 
equation L2.7 is a test to see if the equation is self adjoint or not. 
If it is self adjoint, then no transformations are necessary because 
then the equation can be written directly in the form of LI. 36 . 
Application of the self-adjoint test conditions LI. 34 and LI. 35 
results in equation L2.8 which shows that equation L2.7 is not self 
adjoint and that transformations will be required. It is interesting 



to note at this point that q^ = 0 is satisfied for the simplified 
version of the LPCM solved analytically in Chapter L4 showing that 
it is already self adjoint. 

Since equation L2.7 is non self adjoint, the transformations 
L1.37i LI. 38, and LI. 39 can be applied, as in equations L2.9» L2.10, 
and L2.ll, resulting in the self-adjoint equation LI. 40. With these 
functions defined, the second transformation, to the Liouville Normal- 
Form equation, can be utilized as in equations L2.12 and L2.13 
resulting from equations LI. 44 and L1.45» Application of transformation 



q x - 0 , but q 1 j 0 in PLPBV 



L2.8 



?(x) - exp [(x-a) Iq/Pq + (x 3 - 3 - 3 )^/^] L2 »9 



Q(x) - P(x)q-;/p o 



L2.10 



R(x) * P(x)/p Q 



L2.ll 



X(x) - W(x) ^Tp^ exp [-(x-a)q Q /2p o - 

(x 3 -a 3 )q 1 /4p o ] L2.12 



z = (x-a)/ V~P^ 



L2.13 
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relations LI. 47, LI. 48, and LI, 49 results, after some manipulations, 
in equations L2.14, L2.15, and L2.16. Inverting L2.13 and substitut- 
ing for x in equation L2.14 gives M(z) in L2.17, 



M(x) •* 


V 3 * 


2q 0 q-,x + q 1 2 x 3 ^ 3q T 
4 Pq ~T~ 


L2.14 


c = (b- 


o 

■ 




L2.15 


X = K 3 






L2.16 


M(z) - 




2q o qi ( z + a) + q., 3 


( v iT z + a) 3 






4? o 





3^1 

+ — L2.17 

The next step in the transformation to the Liouville Normal- 
Form equation is the evaluation of the boundary conditions in LI. 50 
through LI. 55* These relationships are given in equations L2.18 
through L2.22, 



h.! ~ ” A 12^ q 0 * q l b >/ 2p o 


L2.18 


D 12 “ A 12/ ^ 


L2.19 


D 21 = A 21 ~ A 22 (q o + q l a )/ 2 Po 


L2.20 


D 22 = A 22^ ^ 


L2.21 


D 11 I) 22 " D 12 B 21 = ^ A 11 A 22 " *21*12^ ^o 




3 

- A 12 A 22 q 1 (b-a)/2p/ / 0 


L2.22 
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The final transformation which must be made is the conversion 

of the Liouville Normal-Form equation to the Fredholm II-Integral 

Equation following the steps in Appendix Al. However, all of the 

functions and constants in equations A1.20 through A1.24 are now 

known in terms of PLPBV, so the transformation is complete. With 

these known transformation relationships and with the solution 

technique of Appendix A3, a digital computer subroutine can be 

written for computation of the solutions X n (x), the eigenfunction 

expansion, and the computation and plotting of the final solution 

u (x,t) . 
m 

L2.5 SOLUTION OF THE TRANSIENT PROBLEM; COMPUTATIONAL ANALYSIS 

Once a complete analytical analysis of any specific subcase of 
the LPCM has been accomplished, the next step is to write a computer 
program to evaluate the solutions utilizing analytical transformations 
and integral equation solution techniques. Such a computer program 
will require a significant amount of work to write and computation time 
to check out. However, it is expected that, in the end, the program 
will require far less computation time than currently available programs 
which solve the discrete -plate equations. It is expected that a 
further significant savings in time would result from the design of a 
special purpose computer to solve a model 'of sufficient degree to 
attain desired accuracy for a specific column or type of column. 

All computer programs utilizing the methods of this thesis will, 
by the nature of the method, follow a certain format. This subsection 
presents a description of the computation steps, in flowchart form, 
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which this solution method uses. This flowchart is presented in 
Figure L2.3, and a Fortran IV program listing the details of the block 
steps in Figure L2.3 for the special case of FLPBV is presented in 
Appendices A5 and A 6. 

The PLPBV subroutine presented in Appendix A6 is intended as a 
description of the steps necessary for solving LPCM problems. Sub- 
routine PLPBV was compiled and executed for several very simplified 
cases but would require a great deal of work to be useful in general. 

It is estimated that PLPBV would require between 100 and 200 man-hours 
of programming time and 1 to 2 hours of computation time to perfect, 
to consider all special cases, and to debug. The primary benefit 
from this labor would be the fact that PLPBV only requires about 15 
seconds of computation time, and thus, complete column solutions using 
PLPBV could be generated in less than one minute for 20 spatial points 
and 9 times. 

This chapter has presented a detailed description of the steps 
necessary to practically apply the general integral equation solution 
technique developed in Chapter LI to a special subcase, PLPBV, of the 
LPCM. These steps sire described analytically by applying the trans- 
formation relations to PLPBV and are described numerically by a flow- 
chart of the necessary computation steps and by a very basic Fortran IV 
digital computer subroutine showing some of the details of the computa- 
tion steps. The next chapter describes the steps necessary to determine^ 
the LPCM for a specific distillation column. 
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Figure L2, 3 (Contd., ) 
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Figure L2.3 




- COMPUTATION STSF3 IN SUBROUTINE FLPBV 
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CHAPTER L3 



DETERMINATION OF THE LPCM TOR EXAMPLE DISTILLATION COLUMNS 

The Linear Polynomial -Coefficient Model (LPCM) has been developed 
and solved in this thesis as a continuous-spatial model of the transient 
composition behavior of a binary plate distillation column. The purpose 
of this chapter is to show the steps required to determine the LPCM 
for given distillation columns and to discuss the approximations 
involved in developing the LPCM. Two columns having the same feed and 
output compositions but different reflux rates and numbers of plates 
are used as examples to show the development of the LPCM. This chapter 
represents an equating of Figure M^-,2 to Figure Ll.l in terms of the 
two example columns. 

L3.1 GENERAL STATEMENT AND DISCUSSION OF THE MODELING STEPS 

The first step in the determination of the LPCM for any distilla- 
tion column is to determine the slope of the equilibrium curve, m(x), 
as a continuous function of x. The function m(x) is defined in equation 
M4»5 and is presented here in detail as equation L3»l« Based upon the 



»(x) = 3f(u)i 



a 



u-u^x) 



[1 + (cr-lju^x)^ 



L3.1 



initial steady-state u. (x.) represented by the McCabe-Thiele diagram, 
the function m(x) can be expressed exactly at only a discrete number 
of points. The major approximation involved in determining the LPCM, 
aside from the original linearization of the CSS, is the expressing 
of m(x) as an n-th degree polynomial in x, as in equation L3#2. 
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, v n 4 

m(x; = Z m. x «* m + ra-, x + m 0 x 3 + L3.2 

1=0 1 o -L d. 

The continuous-spatial steady-state model has not been explicitly 
emphasized in this thesis because the direct approximation of the 
discrete steady-state equilibrium curve slope seems to offer a more 
direct approach. The continuous steady-state model is equivalent to 
Figure M2.1 with all time derivatives set equal to zero* If this model 
were solved for the continuous steady-state u^(x), then the continuous 
m(x) could be found from equation L3.1. This is impractical for the 
LPCM determination because the approximations necessary to develop and 
solve the continuous steady-state model introduce more error than the 
direct determination of m(x) from the discrete steady-state model. 

Thus, continuous steady-state models are not used in this thesis and 
are very seldom found in the literature. 

The accuracy of the LPCM is expected to depend directly upon the 
accuracy of m(x). Theoretically m(x) can be determined to any desired 
accuracy using a polynomial in equation L3.2 of sufficiently high 
degree and using either point-by-point or least-squares techniques. 

The integral equation solution technique of Chapter LI is valid for 
any degree polynomial, but as the degree gets higher than one the 
transformations become overwhelmingly complicated and the programming 
time necessary to implement the transformations and consider all of 
the special bases becomes an order of magnitude larger. Surprisingly 
enough, the computation time, once the proper programs axe written, 
is not expected to increase significantly over that for PLPBV in 
Chapter L2 because the basic steps of Figure L2.3 remain nearly the 
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same and all of the special cases are parallel paths for computation. 

It is expected that models of degree one or two will prove to be 
adequate to specify the column behavior sufficiently for the purpose 
of control. 

Once the continuous function m(x) has been determined, the model 
derivation is essentially complete. All that remains is to calculate 
the polynomial coefficients p^ and q^ and the boundary condition 
constants A^, The complete expressions for evaluation of these 
constants in terms of a second degree m(x) are given in Figure L3« 3, 

L3.2 APPLICATION OF THE MODELING STEPS TO TWO EXAMPLE COLUMNS 

The McCabe-Thieie diagram for a five-plate binary distillation 
column is presented in Figure L3»l and the corresponding diagram for 
an eleven-plate column is presented in Figure 12,2 in Section l(l). 

These two diagrams have been designed to have the same input and output 
compositions to simplify the calculations and to provide for easy 
comparison. The discrete values of m(x) have been calculated using 
equation L3.1 and are presented numerically in Table L3.1 and Table L3»2 
and presented graphically in Figure L3«2. 
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2.31 
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0.36 



Table L3,l - THE DISCRSTE m(x) FOR A FIVS PLATS COLUMN 
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Figure L 3.1 

MCCA3S-THIBLE DIAGRAM FOR 
A FIVE PLATE COLUMN 
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Equilibrium Curve Slope as a Function of 
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Model Equation 

lL [?i(x)u] + 5_ [P ? (x)u] = du 



0X 



B x 



at 



Where: 



P l( x ) = P 0 + PjX + p 2 x 2 = [(B +m o ) + n^x + m 2 x 2 ]/2 

Thus: p = (B + m )/Z 

?! = n]/2 
p 2 " m 2/ 2 

p 2^ x ) = To + 1i x + <1 2 X 2 = [(S " m 0 ) - m 1 x - m 2 x a ] 

Thus: q Q = (B - m o ) 

*1 = ~ m l 
q 2 " m 2 

Both P^(x) and P^x) must be expressed for both the upper 
and lower sections of the column using B u and B^. 

Boundary Conditions 

Upper A in = fdm(x) - m(x) + ? 

b-1.0 J ” L |_ dx J 



a=x. 



x » 1.0 



An = m 2 “ m o + 1 
+ m, + m 0 

a-;2-o 

A 2i = “C 3 U ( m i + Pm 2 x f) + F A] 

k 2Z =1 ~ m o‘ m l x f “ m 2 X f 3 

k = 0 

41 - -F/V 



Lower An = — + 2m 2 x f ) + F/v] 



b=x 
a=0^0 



hz = i " m o - Vf - m 2 X f‘ 

Ano “ 0 
A i4 = -F/V 

-I B ""° 

22 

A = p- 1 - 

A 23 ^ 

l-J - 0 

A 24 ° 

Feed Condition F/V = (B- - B )/q 



Figure L3.3 - THE COMPLETE LPCM WITH SECOND DEGREE POLYNOMIALS 
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Table L3.2 - THE DISCRETE m(x) TOR AN ELEVEN PLATE COLUMN 



The next step in the modeling procedure to develop the LFCM for 
these two columns is to determine the coefficients of a polynomial, 
m(x) , of sufficient degree to meet accuracy requirements based upon 
the discrete points given for m(x). In the case of the 5"Pla-te column 
a polynomial of degree 5 could be found which would exactly match each 
of the 6 points given, and for the 11-plate column a similar polynomial 
of degree 10 could be found. However, for both of these columns a 
2nd-degree polynomial is expected to be accurate enough. Two second 
degree polynomials have been found as examples using the points, 
x = (0.0, 0.2, 1.0) , for the 11 -plate column and the points, 
x = (0,0, 0,6, 1.0), for the 5 -plate column. The resulting m(x) 
polynomials are given by equation L3*3 for the ^-'pla.ie column and 
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L3.4 for the 11 -plate column. 

m(x) = 2.31 - 4.1?x + 2,22x 3 5-plate L3.3 

m(x) = 2.31 - 5*32x + 3»37x 2 11 -plate L3.4 

Once the function m(x) has been determined to the desired accuracy, 
the complete LPCM for the column can be calculated from the equations 
of Figure L3»3» The 2nd-degree m(x) polynomials in equations L3.3 end 
L3.4 result in the complete LPCM for the 5—plete column in Figure L3»4 
and the complete LPCM for the 11-plate column in Figure L3.5* It should 
be emphasized that the approximations and calculations used to obtain 
equations L3»3 end L3.4 were chosen to greatly simplify the calcula- 
tions. Normally, the polynomial coefficients for m(x) should be deter- 
mined specifically for the region of application of the model equation 
instead of using an overall m(x) as was done in equations L3.3 and 
L3.4, The column may also be modeled using more sections than the 
upper and lower sections used in this thesis. 

L3.3 GENERAL COMMENTS AND SUMMARY 

The real utility of the LPCM results from the fact that as the 
number of plates in the column increases, the model complexity remains 
nearly constant for realistic approximations of m(x). The computation 
time and complexity of discrete models, however, increases rapidly as 
the number of plates increases because the dimensions of all of the 
matrices in the model increase with the number of plates. The LPCM is 
valid for components with non-constant relative volatilities as well 
as for constant relative volatilities, but only the comstant relative 
volatility equilibrium curve has been used in this thesis. Also, the 
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5-PLATE COLUMN 




Model Equation 
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Figure L3.4 - THE COMPLETE LPCM FOR 


THE 5-PLATE COLUMN OF FIGURE L3.1 


1 



lpcm can be used for multiple feed and sidestreams and can be applied 
in any manner desired to smaller sections of the column* Each sub- 
section of the column can then be modeled using the LPCM. 
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11 -PLATE COLUMN 



Model Equation 



C p i( x ) u l + Z p 2 ( x ) u l = <Ll 

p a* 2 at 



Where: 



i: P (x) = p + p x + p x s 

*i(x) = q° + q* x + q^x 3 

Upper Section 0.5 £ x £ 1.0 Lower Section 0.0 £ x £ 0.5 
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Figure L3. 5 - THE COMPLETE LPCM FOR THE 11 -PLATE COLUMN OF FIGURE 12,2 



The two models in Figures L3.4 and L3.5 are not solved in this 
thesis. Their solution will require the development of a computer 
program. based upon the integral equation solution technique of Chapter LI 
and using the computation steps of Chapter L2. This will require an 
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extensive amount of analytical analysis, computer programming, and 
testing. The result of this labor is expected to be a general, very 
fast, solution program for the transient behavior of a large class 
of distillation columns which describes the column behavior accurately 
enough for control applications. Extensive research would be required 
to justify this, however. 

This chapter presents the steps necessary to derive the LPCM for 
general distillation columns. The details of these steps are then 
applied to a 5~Pla-te and an 11 -plate column and the resulting LPCM's 
presented. The next chapter solves analytically and computes numeri- 
cally the solutions to some simplified LPCM’s in which m(x) = m Q , a 
constant. 
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CHAPTER L4 



ANALYTICAL SOLUTIONS TO APPROXIMATED COLUMN EQUATIONS 

The purpose of this chapter is to analytically solve and numer- 
ically calculate the total solutions to four simple subcases of the 
Linear Polynomial-Coefficient Model (LPCM) defined in Figure Ll.l, 

Ihe reason for doing this is to provide some examples of simpler 
solutions upon which to base an understanding of solutions to more 
sophisticated versions of the LPCM. The procedure to be followed 
will be to define a simplified version of the LPCM, state the 
problems to be solved, present the analytical solutions to the 
problems, present a computer program used to sum the series, and, 
finally, to present graphically the simplified model solutions. 

L4.1 A SIMPLIFIED LPCM 

A greatly simplified LPCM is presented in Figure L4.1 by 
equations L4.1 through L4.5. Ibis model represents one step below 
the more complicated model presented in Chapter L2, This simplified 
LPCM appears in many places in the literature, only it is known by 
the different names and uses listed in Table L4.1, for a representative 
sample of the literature. This model is usually solved by analytical 
methods in the literature and used as an example of the application 
of the method of separation of variables. 

L4.2 APPROXIMATE EQUATIONS TO BE SOLVED 

Tnere are four specific problems involving the simplified LPCM 
which will be solved analytically in this chapter. These problems 
are presented in Figure L4.2 and Table L4.2, Two of the problems 
involve the solution for the top step response of the heat equation 
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for the two cases of zero bottoms composition and zero bottoms outflow. 
The other two problems solved are the Taylor Diffusion Model (TDM) 
equation for positive and negative values of the center constant 
A visualization of what is happening in terms of distillation 
column concentration or, by analogy, temperature in a rod is presented 
in Figure L4.1. If the response desired is the total response to a 
step change in feed composition, then the equations solved in this 
chapter apply directly to the bottom half of the column in the range 
0 < x £ 1, where x = 1 is the location of the feed tray. The heat 
analogy would be to visualize a rod with length in the range 0 < x < 1 
with a step change in temperature at the top. 
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L4.2 


P 2 (x) . 4. 


L4.3 
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L4.5 


Figure L4.1 - A SIMPLIFIED LPCM 
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AUTHOR (S) REFERENCE EQUATION NAME OR USE 



Crank 


C-8 




Diffusion Equation 


Gould 


G-3 




Taylor Diffusion Model (TDM) 


Jackson and Pigford 


J-l 




Transient-Diffusion Equation 


Lee 


L-ll 






Pollock, Brown 
and Dempsey 


P-5 


> 


Convective Transport Equation 


Stone and Brian 


S-13j 






Bennett and Meyers 








Boyce and DiPrima 


B-8 






Berg and McGregor 


B-21 






Jury 


J-2 






Powers 


P-7 


> 


Heat Equation 


Sagan 


S-10 






Tsang 


T-5 






Woodle 


H-l^ 







Table L4.1 - LITERATURE NOMENCLATURE AND USAGE 
OF THE SIMPLIFIED LPCM 
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Column or Rod 
Visual 




Column Concentration Chart 
Transient Distributions 
Graphical 



Figure L4.2 - VISUAL AND GRAPHICAL SOLUTION FORMAT 



Equation Boundary Conditions 



Name 
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Heat 
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1 


0 


1 
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UA2 


a a u * u, 

XX t 


Heat 


a 3 


0 


1 


0 


0 


1 


1 


UB1 


a a u - u « u ± 

XX x t 


TDM 


tt 3 


-1 


1 


0 


1 


0 


1 


UC1 


nr a u + u *» u, 
U XX x t 
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a 3 


+1 


1 


0 
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Table L4.2 - L?CX PROBLEMS TO 3E SOLVFJ) ANALYTICALLY 
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L4.3 ANALYTICAL SOLUTIONS TO THE SIMPLIFIED LPCM EXAMPLES 



This subsection develops the infinite series representations 
of the analytical solutions to problems UA1, UA2, UB1, and UCl which 
have been summarized in Table L4.2. The manipulations leading to 
these solutions represent no more than a series of mathematical 
exercises. The details of these solutions are included in this 
thesis for the sake of completeness and for the convenience of a 
reader who may want to relate the solution techniques used for the 
general LPCM solution in Chapter LI to the much more familiar 
techniques used in this chapter. 

L4.3.1 SOLUTION OF UA1 

The total solution to UA1 is given by equation L4.1j the 
development is presented in the steps below. 

Statement of Problem t Step L4.1 



x 



A u c U-i (t) 




1 



u(l»t) - U^t) 
u(o,t) = 0 
u(x,o) * 0 



0 



u=0 



■>*t 



Transient plus Steady State: 



Step L4.2 



u(x,t) * u s (x) + u T (x,t) 
a2u sxx " 0 aSu Txx “ u Tt 

u g (0) * 0 u T (o,t) » 0 u T (x,o) - -u s (x) 

u s (l) = 1 UpCl.t) * 0 Uip(x t oo) = 0 



Steady State Solution: 



Step L4.3 



u s (x) - x 
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Transient Solution: 



Step L4.4 



u T (x,t) - X(x) • T(t) 

T n (K) - exp (-K a t) 

X n ( x ) - A^cos K x + Ag sin K x 

a a 

u<p(o,t) - 0 gives A^*> 0 

u^(l,t) - 0 defines eigenvalues sin K - 0 

a 



gives K = n-na 
n 

Eigenfunctions are: 
ajj sin (niTx) exp (-K n 3 t) 

Expanding to meet the initial condition: 

OO 

u^x.O) - -u (x) = -x = E a sin (nrrx) 
j. c> n^o n 

Using the orthogonality property of eigenfunctions: 
(-x) sin (rorrx)dx = E a_ sin (m^x) sin (n-n-x)dx 



HUT 2 



Eigenfunction constants : 
.n 



a n - 

n<rr 

Transient Solution: 

OO 

u T (x,t) » 2 2 (-1) sln(nrrx) exp [-(n<na) 3 t] 

** ott 

Total Solution to UA1: 




L4.1 
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L4.3.2 SOLUTION OF UA2 



The total solution to UA2 is given by equation L4.2, the 
development is presented, below. 

Statement of Problem i 



<* u xx * H = 0 

u(l,t) - u. x (t) 

u x (o,t) - 0 



u « (t) 



•u - 0 



u x “ 0 



a sxx 

“.( 1 ) * 1 



Step L4.1 



t 



u(x,o) “ 0 

Transient plus Steady State s Step L4.2 

u(x,t) - u a (x) + u^x.t) 

0 a3u Txx = u Tt 

Uip(l.t) = 0 Urp(x. o) - -u s (x) 

u sx (0) = 0 u Tx (o,t) = 0 u T (x,oo) = 0 



Steady State Solution 
u s (x) = 1 
Transient Solutions 



u T (x,t) 


« X(x) • T(t) 




T »w - 


exp (**K 2 t) 




- 


cos K x + £ 


sin X x 


a 


CL 


<1^(0, t) 


= 0 gives A£ 


- 0 


U T fl,tJ 


= 0 defines eigenvalues 




gives K n 


=* nna 



Step L4.3 
Step L4.4 



cos K = 0 
a 



Eigenfunctions are: 

Oyj cos (mrx) exp (-K n 3 t) 

2 

Expanding to meet the initial conditions 
u T (x,0) - -u s (x) = -1 = E =Q a n cos (ngx) 
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Using the orthogonality property of these eigenfunctions! 

OO 1 

/ o (-1) cos (rorrx) dx - Z a„ cos ( onx) . 



cos ( nirx) dx=*> sin (rnv/2) =-a_ 

2 rim/ 2 £ 

Eigenfunction constants: 

a « -2 sin (nm/2) 
n irm/2 

Transient Solutions 

OO 

u„(x,t) = -2 E sin (nn/2) cos (nmx/2) exp ^-(nna/2) 3 t3 
1 n “ 1 



Total Solution to UA2: 



OO 



u(x* t) - 1 - 2 Z 
n»l 



sin(nTr/2) cos (nTrx/2) exp f-(nna/2) a t] 
ntr/2 

L4.2 



L4. 3. 3 SOLUTIONS OF UB1 and UC1 

The total solutions to problems UB1 and UC1 are given in 
equations L4. 3 and L4.4, respectivexy. The development of the UB1 
solution is presented below; the UC1 solution development is identical 
to UB1 except for the obvious sign change. 



Statement of Problem UBlt 

a 2u xx " u x “ u t 
u(l,t) » U_ x (t) 

u(o,t) = 0 
u(x,o) *» 0 



Step L4.1 



u 



U.i(t) 



u 



0 



-3* t 
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Transient plus Steady States 



Step L4.2 



u(x,t) = u s (x) + u^x.t) 

aSu sxx “ u x " 0 a3u Txx * “TX “ “Tt 
u s (°) = 0 u T (o,t) = 0 u T (x,o) 

u g (l) = 0 u T (l,t) =* 0 u T (x,<») 

Steady State Solutions 

u s (x) » A + B exp (x/a 3 ) 

u g (o) » 0 gives A + B = 0 

u g (l) =* 1 gives A + B exp (l/a 3 ) » 1 

Solving the two A, B equations gives 

« (x) = 1 “ exp (x/g 3 ) 

1 - ex P (1/a 8 ) 

Transient Solutions 

u,p(x,t) » X(x) • T(t) 

T n (t) - exp (-K 2 t) 

X n (x) = (A-jCos Dx + A^ sin Dx) exp (x/2a 3 ) 

u T (o,t) « 0 gives A^ » 0 

u^(l,t) » 0 defines eigenvalues sin D =» 

where s D„ 2 = Jta 3 !^ 3 - 1 = (n-n) 3 
2a 2 

solving for the eigenvalues 

K 3 =■ 1 + (nrtt) 3 

4a 3 

Eigenfunctions are: 

a^ sin (mfx) exp (x/2a 3 ) exp (-K n 3 t) 

Expanding to meet the initial conditions 

u T (x„o) - -u g (x) » - [1 - exp (x/a 3 )]/C 

= £ „a_ sin (mix) exp (x/2a 3 ) 
n=l n 

where: C » 1 - exp (l/a 3 ) 

-153- 



- -U g (x) 

- 0 

Step L4.3 



Step L4.4 



0 gives D 

n 



Using the orthogonality with weighting function property 



of the eigenfunctions s 

S Q £exp (x/a a ) - 1^ exp (-x/2q a ) sin(mTTx) dx 

c 

= E a sin (ruTTx) sin (nrrx) dx 
n=l n ° 

(-1) (mrr) exp (-l/2a a ) = a m 

-f- 

Eigenfunction constants: 

[M 3 + (i/2a a ) 3 ] 

Transient Solution: 

u T (x,t) = 2 £ r(x-; 1 )/2a a 1 sin (nnx) 

n L(nrr) a + (l/2a 3 ) 3 ] 

exp [-(l/^a s + a 3 n 3 ™ 2 ) t] 

Total Solution to UBli 



u(x,t) - 


"l - exp (x/a a )~i + 2 exp T-t/^a 3 + (x-l)/2a a l» 




[l'- exp h/a4j 


E. (-1) (nrr) exp f- (anTr) a tl sin (nrrx) 

n ^LT^FVTl72j)h 




L4.3 



Total Solution to UC1: 
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L4. 4 NUMERICAL CALCULATION AND GRAPHICAL PLOTS OF THE SOLUTIONS 



Each of the solutions to the four problems of Table L4.2 contains 
an infinite series in eigenvalues and eigenfunctions. This section 
describes briefly a digital computer program which was written and 
utilized to numerically calculate the values of u(x,t) for a “ 1.18 
at ten spatial points and nine time points. These values were 
determined for the purpose of making graphical plots of the solutions. 

The computer printouts and description of the program are given 
in Appendix A4. The computer output plots have been used to make the 
following graphical presentations of the four solutions which are 
presented in Figures L4.3 through L4. 6. It is helpful to visualize 
these solution curves in terms of Figure L4. 2. Each solution curve 
from left to right represents a spatial distribution at a given time. 

The curves can be seen to approach their steady state values as, for 
example, the linear distribution in UA1 or the curved exponential 
distribution in UC1. 

The digital computer program in Appendix A4 computes all four 
of the solutions to three decimal place accuracy. It was found that 
at all spatial points and time points the number of eigenvalues and 
eigenfunctions, i.e. terms in the series, necessary to achieve three 
decimal place accuracy was in all cases less than ten (10). 

L4. 5 SIMILARITIES BETWEEN CHAPTER LI SOLUTIONS AND CHAPTER L4 SOLUTIONS 

It is interesting at this point to investigate the similarity 
between the standard analytical solution methods of this chapter and 
the general solution technique presented for the LPCM in Chapter LI. 

This relationship can best be seen by relating each Step in this 
chapter to the corresponding equations in Chapter LI. 
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Ihe problem statement Step L4.1 or, more generally, Figure L4.1 
corresponds to Figure Ll.l, the statement of the LPCM. The transient 
and steady state conditions are similar to equation LI. 15 and the 
boundary conditions LI. 8 through LI. 13. In both cases, the steady 
state solution, Step L4. 3 or LI, 5, LI. 8, and L1.9» represents the 
solution to a two point boundary value problem. This is easily 
found in the simplified equations but requires a computer program 
(such as LBVP) for the general case. 

The transient solution Step L4.4, is, of course, the major step 
in solving the simplified examples as it is, also, in solving the LPCM 
in general. In the simplified examples, the eigenfunctions X n (x) 
are expressible as closed-form functions with very obvious properties 
which can be used to satisfy the boundary conditions and to immediately 
solve for the eigenvalues K n . In the general LPCM, however, once the 
equation LI. 25 has spatially varying coefficients, the eigenfunctions, 
in almost all cases, cannot be expressed in closed form but must be 
represented by infinite series. The infinite series representations 
of the eigenfunctions can be found by a number of techniques, such as 
the power series Method of Frobenius (H-12, Ch. 4), but these methods 
have several major disadvantages in this application. 

The first major disadvantage to infinite-series eigenfunction 
representations is that their properties in terras of satisfying 
boundary conditions are usually not obvious. Secondly, these methods 
lend themselves best to specific example solutions and become very 
difficult to talk about in terms of a general analytical solution 
technique because of the “special cases" involved. Finally, the 
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eigenfunction expansion to meet the initial condition is analytically 
difficult, in general, when these infinite-series eigenfunctions are 
involved. For these reasons, two other analytical techniques for 
finding the eigenvalues and eigenfunctions have been investigated, 

These two techniques are the Prufer Substitution (B-14, Ch. ll) 
and Fredholm II - Integral Equation theory (Appendix Al). The Prufer 
Substitution technique is a transformation applied to the second- 
order differential equation L1.25 which results in two first order 
ordinary differential equations whose solutions define the eigenvalues 
and eigenfunctions. The resulting first order equations are very 
useful in showing the properties, such as existence, ordlnality, 
separation, orthogonality, etc,, of the eigenvalues and eigenfunctions, 
but the differential equations are analytically un tractable in terms 
of specific solutions and did not seem to lend themselves to consistent 
approximations. For these reasons the Prufer Substitution has been 
discarded as a possible analytical, solution technique and the integral, 
equation technique has been developed in this thesis. 

The final part of the transient solution Step L4.4 is the eigen- 
function expansion to meet the initial condition L1.12. The key to 
success in this expansion for the simplified cases is the orthogonality 
of the eigenfunctions with respect to a weighting function over the 
same interval as the boundary conditions are applicable. This 
orthogonality condition is expressed by equation L4.5. Of course, the 
sine and cosine functions of the simplified LPCM examples satisfy L4.5. 
The real utility of the eigenfunction orthogonality condition in the 
simplified examples lies in the fact that the constants A fi can be each 

explicitly expressed in terms of n, 
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n 



4 “ (x > \(*)X„(x)d* - 6 m C n L4.5 

where: ® constants / 0 

6 “0 for m ^ n 

mn ' 

6 " 1 for m =» n 

mn 

One may well wonder why the orthogonality condition was not 
utilized in solving for the constants C in equation LI. 57, when this 
is one of the most important properties of these functions. The reason 
is that the functions X R (x) are the result of transformations LI. 44 
and LI, 45 applied to the solutions W^(z) obtained from Appendix A2 
or Appendix A3 and, as such, are not easily expressed in general 
form, let alone integrated. For specific cases of the LPCM, if the 
functions X^(x) are expressible in a form easily evaluated by inte- 
gration, then the set of equations L1.59» only with the orthogonality 
constants instead of the functions evaluated at points x., will be 
individually solvable because the matrix X will be diagonal. In that 
case, a better solution technique would be to use the orthogonality 
condition rather than the discrete-point approximate solution technique 
employed in equation LI, 58 . 

The fact that the functions X (x) satisfy the orthogonality 

n 

condition, in the general case, is a result of the characteristics 
of the original eigenvalue problem LI. 25* If the transformation 
relations used in subsections LI. 5 and LI. 6 satisfy continuity con- 
ditions and result in the regular Sturm Liouville problem LI. 40, then 
the kernel (A1.17) of the integral equation is guaranteed to be con- 
tinuous, In addition, it has been shown by Lovitt (L-8, pp. 181-182) 
that for the special case of R(x) =* 1 in equation LI. 40, the kernel 
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A1.17 is symmetric and satisfies equation L4.6. The kernel symmetry 



K(z,s) - K(s,z) L 4,6 

guarantees the existence and reality of the eigenvalues and the 
orthogonality of the eigenfunctions. Of course, in the simplified 
LPCM examples of this chapter the function R(x) io always one (l). 
In summary, this chapter presents analytical and graphical 
solutions to four simplified cases of the LPCM. The techniques for 
finding these solutions is then related to the general, integral - 
equation solution technique for the LPCM presented in Chapter LI. 
The next Chapter L5 presents a brief discussion of several alterna- 
tive solution techniques which might be applied to the LPCM. 
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CHAPTER L5 



OTHER SOLUTION TECHNIQUES APPLICABLE TO THE LPCM 

This thesis presents the LPCM as a dynamic model of the composi- 
tion behavior of a binary plate distillation column and suggests an 
integral equation solution technique for it. The speed and precision 
of the integral equation solution technique has not been demonstrated 
completely in this thesis, however. Therefore, this chapter describes 
briefly the basic principles behind several other numerical methods 
for the solution of parabolic partial differential equations which 
might be applied to the LPCM for the purpose of comparing numerical 
solution accuracy and speed with the integral equation method. 

Numerical methods for the solution of ordinary and partial 
differential equations are extremely useful in the study of distilla- 
tion columns since both the discrete-plate and the continuous-spatial 
equations are impossible to solve analytically except in the simplest 
cases. Even when analytic methods are applicable, it is often the 
case that the solution is expressed in the form of an infinite series 
rather than in closed form. When this is the case, it may sometimes 
be less time consuming to apply a numerical method directly than to 
evaluate a series to some desired degree of accuracy at each point 
of interest. This was not true in the simplified cases of Chapter L4 
because each of the series converged so rapidly, i.e. less than ten 
terms. 

There are a large number of numerical methods which can be applied 
to any particular distillation column model, such as the LPCM. The 
choice of a particular one may depend on the form of the LPCM for any 
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particular column. This chapter presents the basic principles behind 
the explicit numerical methods used for solving parabolic partial 
differential equations, applies these methods to the LPCM, and presents 
the resulting recursion relations, 

L5.1 NUMERICAL SOLUTION TECHNIQUES FOR PARABOLIC PARTIAL DIFFERENTIAL 
EQUATIONS 

Parabolic partial differential equations are usually solved by 
the use of step-by-step finite difference methods (B2.4 - Numerical 
Solution Techniques) . This method consists essentially of defining 
a regular (usually rectangular) mesh and replacing the differential 
equation by a difference equation defined on the nodes of the mesh. 

Most of the variations in types of methods result from different 
derivations of the difference equation from the differential equation. 

In the distillation column models described in Section 2 (m), 
the solution sought is a function u(x,t) where both u and x are in 
the 0 to 1 interval and t is greater than zero. A possible mesh for 
use in equations of this type is presented in Figure L5.1. In this 
mesh the solution domain is divided into intervals of width h in x and 
1 in t. Ihe two increments h and 1 need not be the same and in some 
cases may be changed over different portions of the domains. The 
ratio of these two increments is constrained by stability requirements 
for any particular equation and method of solution. 

Once the mesh is set up the next step is to consider u(x,t) only 
at the mesh junction points as in equation L5.1. The parabolic LPCM 

u(x,t) “ u(ih, jl) = u, . L5.1 

t 
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can be solved by either explicit or implicit methods. Explicit methods 
use only previous and present values of u, such as u. u, • , , • • •, 
to compute the future value u i } j+i • Implicit methods use previous, 
present, and future values of u to approximate the derivatives at the 
present point (ih, jl) which are then used to refine the present value 
of u. Both methods are applicable to the LPCM, but only explicit 
methods will be discussed here. It is expected that explicit methods 
would be adequate for solving most cases of the LPCM, but if not, 
then implicit methods (B2.4 - Numerical Solution Techniques) could be 
employed. Explicit methods are generally easy to implement and rapid 
to execute but are relatively unstable compared to implicit methods. 
Most finite difference approximations are derived from a Taylor 
Series expansion of u about the value u^j. Only as many terms in 
the series are retained as necessary for an accurate solution. A 

Taylor Series expansion in t with only the first order terms retained 
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is presented in equation L5.2. Similarly, Taylor Series expansions 
in x with up to the second order terms retained are presented in 
equations L5.3 and L5.4. 



u . , « u. . + du 
1,3+1 1,3 p 



L5.2 



i ,3 



Vi. j * “i.j + 



h + 3fu 



i* j 



h s /2 L5.3 



i. 3 



if j 



h + 2 s u 

a? 



h a /2 L5.4 



i. j 



L5.2 APPLICATION OF EXPLICIT METHODS TO THE LPCM 

The basic idea is to use stepwise expressions for u and its 
derivatives in the LPCM equation, here presented in expanded form is 
L5.5. In that case the above expansions must be solved for the 



P (x) 4 3 u + P~(x) 3u + P^(x)u = 0u L3.5 

1 a * 3 a* d t 

approximated derivates of u in terms of the previous and present 
values of u. These solutions will depend upon how many terms of the 
expansion are used. For the above expansions the approximated values 
of u and its derivatives are given by equations L5.6 through L5.9. 



u = u 



if j 



3u 

3 X 



i» j 



u ifj : u i-l«j 



L5.6 

L5.7 



3 a u 

3x 3 



i.j 



U i-H , j ~ 2u i.i + u i-l,j 
h 3 



L5.8 
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Substituting these approximate relations into the LPCM equation and 
then solving for u i+1> j gives the recursion relation L5.10 for the 
LPCM, where the polynomial coefficients are evaluated at x = ih. 



U i+1 » j = ^ " A i ~ B i + C i^ u i, j " A i u i-l,j 

L5.10 

" C i U i,j-l 

Where t = 1 - Pyi/P^ 

B. = P 4 h 2 /P 1 

- h 3 /(P 1 l) 

P i (x) are defined in equations LI. 30 
through LI. 33 

Using L5.10 and the previously derived equations a numerical iteration 
scheme can be set up for solving the LPCM along the mesh. 

L5.3 BRIEF SUMMARY OF SECTION 3(L) 

Section 3(L) concentrates on the Linear Polynomial -Coefficient 
Model (LPCM) which is derived in Section 2 (m) to represent the 
transient behavior of a binary plate distillation column. Chapter LI 
presents an integral equation procedure for analytical solution of the 
LPCM in general using Appendices Al, A2, and A3* Chapter L2 applies 
the solution technique to a reduced form of the LPCM having one spatially 
varying coefficient and one constant coefficient, A suggested scheme 
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for a computer program PLPBV to evaluate the solution is presented 
using Appendices A5 and A6, Chapter L3 evaluates the LPCM for several 
cases from the steady-state representation of an example column. 
Chapter L4 solves analytically and evaluates numerically using Appen- 
dix A4, four simplified examples of the LPCM: two cases of the best 

equation and two cases of the Taylor Diffusion Model. The analytical 
solution steps for the simplified examples are then related to the 
analytical steps of the integral equation technique of Chapter LI. 
Chapter L5 discusses briefly several possible numerical methods for 
solving the LPCM and derives the detailed recursion relations to 
solve the LPCM by an explicit numerical method. 
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SECTION 4 



SUMMARY AND CONCLUSION ( S ) 

51 SUMMARY AND CONCLUSIONS 

52 AREAS FOR FURTHER STUDY 



" THE PURPOSE OF COMPUTING IS INSIGHT, NOT NUMBERS" - R. W. HAMMING 

(H-19) 



THE CENTRAL PURPOSE OF THIS THESIS IS THE DEVELOPMENT, PRESENTA- 
TION, SUGGESTED SOLUTION TECHNIQUE, AND EVALUATION OF THE LINEAR 
POLYNOMIAL-COEFFICIENT MODEL (LPCM) OF THE DYNAMIC BEHAVIOR OF A 
BINARY PLATE DISTILLATION COLUMN. THIS SECTION SUMMARIZES THE MAJOR 
ARGUMENTS PRESENTED IN THIS THESIS RELATING TO THE LPCM AND EMPHASIZES 
A LARGE NUMBER OF AREAS FOR FURTHER STUDY USING THE LPCM AND THE 
INTEGRAL EQUATION SOLUTION TECHNIQUE PRESENTED IN THIS THESIS. 



P 
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CHAPTER SI. 



SUMMARY AND CONCLUSIONS 

Ihis thesis is summarized, in the format of the Chapter and 
Appendix Relationship Diagram presented in Chapter C2. 

Sl.l THE MAIN READING PATH OF CHAPTER C2 

Starting from the basic principles of distillation this thesis 
develops a discrete-plate dynamic model for the composition behavior 
of a binary plate distillation column. Ihe fact that discrete models 
characteristically have large solution times leads to the search for 
a faster model and to the investigation of a continuous-spatial 
dynamic model derived by treating the plate number in the discrete 
model as a continuous variable. An investigation of possible solution 
techniques for the continuous-spatial model leads to the conclusion 
that for any hope of an analytical solution, the nonlinear continuous- 
spatial model must be linearized. The representation of the spatial 
coefficients of the linearized continuous model by general n-th degree 
polynomials defines the Linear Polynomial -Coefficient Model (LPCM). 

The basic equation of the LPCM is presented here as equation Sl.l. 

# 2 [P. (x)u"j + 3_ [Po(x)u] « 5_U Sl.l 

a* 1 1 a* a* 

Where i P. (x) and P 2 (x) are polynomials in x. 

Two-point boundary conditions and step initial 
conditions are specified for u and u x# 

Ihe steady-state solution u (x) of the LPCM is obtained by solving the 

s 

two-point boundary value problem. The transient solution is derived 
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by converting the spatial ordinary differential eigenvalue problem 
resulting from separation of variables to a Liouville Normal -Form 
equation and further transforming to a homogeneous Fredholm II integral 
equation. The solutions of the integral equation then define the 
eigenvalues K i and the eigenfunctions X i (x). An eigenfunction expan- 
sion to meet the initial condition then defines the eigenfunction 
constants C^, The total solution to the LPCM is then given by equation 
SI. 2. 

u(x,t) « u (x) + E C.X (x) exp(-K i a t) SI. 2 

s i=l 1 1 1 

The integral equation solution technique used to develop equa- 
tion SI. 2 is described in detail by applying it to a first-degree 
polynomial model (PLPBV) and by suggesting the type of computer 
program which must be used to evaluate the solution. The individual 
steps in the solution technique are described analytically by applying 
the transformation relations to FLPBV and are described numerically 
by a flowchart of the necessary computation steps and by a very basic 
Fortran IV digital computer subroutine showing some of the details of 
the computation steps. Preliminary tests with the PLPBV program 
suggest that it may be possible to generate complete column solutions 
in less than one minute for 20 spatial points and 9 time points. 

The complete model of a binary plate distillation column using 
the LPCM consists of two equations of the form of Sl.l and their corres- 
ponding two-point boundary conditions. The steps necessary to derive 
the LPCM for general distillation columns are presented. The details 
of these steps are then applied, for a second degree model, to a 
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5-plate and an 11-plate column and the resulting models are presented 
as examples. 

Four simplified versions of the LPCM in the form of the heat 
equation and the Taylor Diffusion Model are defined as simplified 
versions of the LPCM. These models axe then solved analytically, a 
computer program used to sum the series is developed, and the solutions 
are presented graphically and numerically. The standard analytical 
solution methods for the simplified models are then compared to the 
general solution technique for the LPCM. 

There are a large number of numerical methods which can be 
applied to any particular distillation column model. Several possible 
numerical methods for solving the LPCM are discussed, and the detailed 
recursion relations to solve the LPCM by an explicit numerical method 
are derived. 

Ihe thesis is summarized briefly, and a discussion of a large 
number of areas for further study is presented. 

SI. 2 THE AREAS OF OVERALL PERTINENCE OF CHAPTER C2 

A bibliography containing 352 references of which 202 pertain to 
distillation column dynamics and control is presented. These references 
are then related to four major areas pertinent to the thesis. Some 

1. General Theory of Distillation 

2. Distillation Column Dynamics 

3. Distillation Column Control 

4. Mathematics and Computation 

general notes, descriptions, and comments pertaining to most of the 
references are included. 
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Several opinions on the philosophical aspects of observing 
reality as related to modeling a distillation column and on the desire 
to achieve profit as related to the cost of separation of components 
are presented. 
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CHAPTER S2 



AREAS FOR FURTHER STUDY 

This chapter presents, in the form of suggestions, a list of 
possible topics for further study or areas for further research in 
the order in which they are encountered by following the main reading 
path of Chapter C2, It must be emphasized that the LPCM and the 
integral equation solution technique are suggested by this thesis 
to be valuable tools for modeling distillation columns, but the vali- 
dation of this suggestion car only come from further study, 

S2.1 SECTION M 

1. Using the general discrete equations developed firom considering 
all four (or part) mass balances per plate, develop continuous-spatial 
models and then LPCM models for the general binary plate column, 

2. Investigate other methods of converting discrete models to continu 
ous models which may give continuous models of greater accuracy, 

3. Investigate analog computer solution techniques for the CSE and 
the nonlinear polynomial -class CSE. 

4. Find in the literature or develop transformations applicable to 
the CSE,, 

5. Apply the Fixed Point Theorem to the CSE, 

a. Show that the CSE is a contraction mapping, 

b. Investigate techniques for finding the fixed point for the CSE 

6. Investigate solution techniques for polynomial-class partial dif- 
ferential equations to see if nonlinear approximated CSE's are more 
easily solved and less accurate than the general CSE. 
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7. Investigate the validity and accuracy of models linearized about 
operating points other than the initial steady-state. 

S2.2 SECTION L 

8. Investigate the validity of the LPCM as a distillation column 
model. using programs developed and based upon the integral equation 
solution technique of this thesis. 

a. Complete PLPBV including all special cases. 

b. Develop main modeling program to start from the equilibrium 
curve and the column physical characteristics and to end up 
with complete LPCM and solutions using PLPBV. 

9. Develop the digital computer programs for the second-degree (and 
higher) versions of the LPCM and use them to compute column solutions. 

10. Develop numerical solution program for the LPCM and compare 
solution time with the analytical solution programs. 

11. Extend the step solution capability of the LPCM to approximate 
the response to any arbitrary input function of time. 

12. Investigate the relative magnitudes of the eigenvalues for several 
variations of the LPCM in the interest of seeing how many are really 
necessary to completely characterize column behavior, 

13. Apply separable kernel solution procedure to the integral equa- 
tion resulting from the LPCM and evaluate the accuracy, 

14. Investigate in detail the properties of the function M(z) for 
several variations of the LPCM. 

15. Apply the integral equation solution technique to the two-point 
steady-state boundary-value problem and combine with the transient 
portion. 
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16 . Develop an LPCM modeling technique for the case when the relative 
volatility is not constant. 

17. Develop an LPCM modeling technique for the case when the Murphree 
plate efficiencies vary and are not necessarily unity. 

18. Investigate in detail the symmetry and continuity properties of 
the integral equation kernel K(z,s) for various functions M(z). 

19. Investigate implicit numerical methods for solving the LPCM and 
compare solution times to the analytical methods. 

20. Investigate the possibilities for partially analytical - partially 
numerical solutions to the LPCM (hybrid solutions). 

S2.3 EXTENSIONS 

21. Develop distillation column control schemes particularly suited 
to using the LPCM. 

22. Investigate the possibilities for partially discrete - partially 
continuous models for distillation columns (hybrid models). 

23. Extend the application of the LPCM techniques to packed columns. 

24. Extend the application of the LPCM to continuous models of multi- 
component distillation columns. 

25. Combine spatial LPCM and time LPCM solution techniques to solve 
partial differential two-point boundary-value problems in which the 
partial differential equation is second order in both time and space. 
The result of the application of separation of variables will be a 
spatial LPCM and a time LPCM. 



- 177 - 



SECTION 5 



APPENDICES (A) 

A1 CONVERSION OF A LIOUVILLE NORMAL-FORM EQUATION TO A 
FREDHOLM II - INTEGRAL EQUATION 

A 2 SOLUTION OF A HOMOGENEOUS FREDHOLM II - INTEGRAL EQUATION 
WITH A SEPARABLE KERNEL 

A3 SOLUTION OF A HOMOGENEOUS FREDHOLM II - INTEGRAL EQJATION 
BY CONVERSION TO LINEAR ALGEBRAIC EQUATIONS 

A4 DIGITAL COMPUTER PROGRAM AND OUTPUT FOR EVALUATION OF 
PROBLEMS IN CHAPTER Lk 

A5 SUBROUTINES USEFUL FOR COMPUTING LPCM STEADY-STATE SOLUTIONS 
A6 SUBROUTINE PLPBV 



"MATHEMATICS POSSESSES NOT ONLY TRUTH, BUT SUPREME BEAUTY — 

A BEAUTY COLD AND AUSTERE, LIKE THAT OF SCULPTURE, • • • SUBLIMELY 
PURE, AND CAPABLE OF A STERN PERFECTION SUCH AS ONLY THE GREATEST 
ART CAN SHOW." - BERTRAND RUSSELL 
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APPENDIX A1 



CONVERSION OF A LIOUVILLE NORMAL-PORN EQUATION TO A 
FREDHOLM II - INTEGRAL EQUATION 

This appendix presents the intermediate steps necessary to the 
conversion of the Liouville Normal-Form equation (B-14) and its 
boundary conditions to a Fredholm II - Integral Equation (L-8.H-13). 
The conversion proceeds via integration of equation Al.l and 
application of boundary conditions A1.2, A1.3# and A1.4 to integral 



equation of the form of Al,5« 

d a W - [M(z)-x]w = 0 Al.l 

dz 5 

D^Wtc) + I^gWCc) - 0 A1.2 

D 2 iW(o) + D 22 W(o) =0 A1.3 

d H d 22 “ D 12 D 21 ^ 0 

V(z) ■ /q K(z,s)W(s)ds A1.5 



Integrating Al.l from 0 to Z twice and making a change of 
variable in one of the integrals yields A1.6, Application of the 

W(z) « Si [M(z-s)-x]w(s)ds + W(o)z 

+ W(o) A1.6 

second boundary condition Al#3 to A1.6 then gives Al*7« The 

W(z) - f [M(z-s)-x]w(s)ds + [l-D 21 z]w(o) A1.7 

°22 

remainder of the conversion procedure consists of determining 
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the constant W(o) and placing the resulting equations in integral 
equation form. 



The first step in evaluating W(o) is to integrate equation 
Al.l from Z to C, giving equation A1.8. Evaluation of A1.8 at Z - 0 

W(e) « W(z) + S° [M(z)-x]w(z)dz A1.8 

z 

suid of A1.7 at Z » C gives A1.9 and A1.10, Placing the boundary 

W(c) o 8(o) + -S° o [M(z)-\>(z)dz A1.9 

W(c) = J* C [M(c-s)-x]w(s)ds + [l - ®21 c]w(o) A1.10 
° D 22 

conditions in the form of Al.ll and A1.12, substituting them into 
A1.9 and A1.10, and equating them gives A1.13, an equation solely 
in terms of W(0). 

W(c) = - °12 W(c) Al.ll 

^1 

o D 

W(o) « - _21 W(o) A1.12 

D 22 

,f C rM(c-s)-x]w(s)ds + [l - ^1 c]w(o) 

0 D 

^22 

- - ^2 [- ^21 W(o) + J“ C [M(z)-x]w(z)dz] A1.13 
»ll D 22 

Solving for W(0) gives A1.14 and A1.15. 

u(o) - X /°[(H(o- 3) + “l2 M(s)) - (1+J2)i] 
c ° Dll “ll 

W(s)ds A1.14 
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A1.15 



C o - ^12 21 + ^21 - 1 / 0 

®U^22 D 22 

If equation A1.14 for W(o) is now substituted back into equation 
A1.7, then equation A1.16 is an integral equation in W(z). This 
equation Al,l6, however, is not yet in the desired form of Al,5. 

W(z) - f [M(z-s)-x] W(s)ds + 1_ q- D 21 z) 

° C o D 22 

S° [(M(c-s) + \z M(s) ) - (l+^L2)x] 

° D 11 D 11 

W(s)ds A1.16 

If the kernel function K(z,s) is expressed as the sum in equations 
A1.17, A1.18, A1.19, and Al. 20, then the conversion to the integral 
equation is complete, 

K(z,s) = K 1 (z , s) + K 2 (z,s) 

s <2 S>z A1.17 

( z , s ) = M(z-s) + 1_ (D 2 2“ D 21 z ) 

C 1 

(D n M(c-s)4C 12 M(s)) 

Cl + (^22~^21 z ^ ^11 + ^12^ Al,l8 

C 1 

K 2 (z,s) « 1_ (D 22 -D 21 z ) [D i;l M(c-s)+D 12 M(s) 

C 1 

-(Dh + Di 2 )a] ai.19 

C 1 “ ^l^ 0 © “ D 12 D 21 “ D ll D 22 + ^Ll^l A 1 * 20 
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Thus, the complete Integral, equation representation of the 
Liouville Normal-Form equation is summarized in equations A1.20 



through A1.24. 

W(z) = J“° K(z,s)w(s)ds A1.21 

K(z,s) » K-^z.s) + K 2 (z,s) A1.22 

K^(z,s) « M(z-s) - x + K 2 (z,s) s < z Al .23 

K 2 (z,s) a 1_ (D22-D21 21 ) C D 11 M ( C " S )+ D 12 M ^ S ^ 

C 1 

- ( D ii +D i2)^ s > z A1.24 



Now, the kernel given in A1.22 must have two properties in order 
for solution procedures to be applied. It must be continuous at s » z 
and it must be symmetric. The continuity aspect implies that Al.23 
can now be written as A1.25 and A1.26, In this case a new function 

K^z.s) = M(z-s) - M(0) + K 2 (z,s) A1.25 

K^(z,s) = \K-j(z,s) A1.26 

K^(z,s) has been introduced, where K 2 (z,s) is given by A1.27 and 
A1.28. Utilizing once again the continuity requirement and the 

K 2 (z,s) => -M(z,s) + M(0) + \K^(z,s) A1.27 

K 2 (z,s) = \K^(z,s) A1.28 

definition of K^(z,s) given in A1.28, then JC^(z,s) can be written 
in terms of K^(z,s) as in equation A1.29. The remaining kernel 
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IfyXz.s) - 1 - M(z-s) + Ko(z,s) A1.29 

W 

Inhere: M(0) / 0 

function K^(z,s) can then be expressed as in equation A1.30, and 
a general test for symmetry can be given by equation A1.31. 

Kt(z,s) = M(z t s) -1+1 - D 0 , z) • 

J m(o) ^ 21 

[D u M(c-s) + D^MCs) - (D^ + D^)] A1.30 
M(0) 

C x [m(zts) - M(0)] + (D 22 - D 21 z) [^(c-s) + D 12 M(s) 

- M(0) (D n + D 12 )] = (D 22 - D 21 s) [D n M(c-z) 

+ D 12 M(z) - (D n + D 12 ) M(0)] A1.31 

The summarized equations for the integral equation representation 
are presented in Figure Al.l. 

W(z) - \ K(z,s)w(s)ds 
K(z,s) - K^(z,s) + K 4 (z,s) 
s < z s > z 

K/,(z t s) » 1 - M(z,s )_ + K~(z,s) s > z 

mToT J 

Ko(z,s) » M(z-s) -1+1 (D ?? - D z) • 

J K(0) C-l 22 21 

^-^(c-s) + D^s) - + D 1? ) n s < z 

*- m(o5 ' 

Figure Al.l - THE FREDHOLM II - HOMOGENEOUS INTEGRAL 

EQUATION 
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APPENDIX A2 



SOLUTION OF A HOMOGENEOUS FREDHOLM II- INTEGRAL EQUATION WITH 
A SEPARABLE KERNEL 

This appendix presents the developments necessary for finding 
the solution to a separable kernel integral equation. The integral 
equation A2.1 and separable kernel A2.2 are combined and a matrix 
solution technique is developed. The technique to be presented is 
standard throughout the literature when it is recognized that the 
following terms are equivalent to "separable kernel" A2.2. 



1. Separable Kernel (H-13) (L-8) 

2. Degenerate Kernel (P-4) (G-?) 

3. Kernel of Finite Rank (S-ll) 

4. Pincherle-Goursat Kernels (T-3) 

5. Riesz-Schauder Equations (M-6) 

W(z) - X / C K(z,s)w(s)ds A2.1 

o 

K(z,s) s J E l a i (z)b, (s) A2.2 



If A2.2 is substituted into A2.1 and the integral and summation 
signs interchanged, then equation A2.3 results. It can be seen that 
the integrated terms depend only upon S and are therefore constants 
A2.4 after integration. Thus, the solution to A2.1 can be expressed 



in terms of the constants as equation A2.5, and the problem now 
becomes that of finding the C . . 



W(z) - A i | l a i (z) b i (s)w(s)ds 

C i “ fo b i( s ) w ( s ) ds 



A2.3 



A2.4 



i* 
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m 

W(z) =■ 



A2.5 



The first step in finding the is the substitution of A2.5 
back into A2.4 giving A2.6. The terms involving the integral of 

a^(s) and b^(s) are then recognized as constants A^.of equation A 2.7, 

3 

which can be calculated from the given kernel function A2.2, The 
resulting equation A2.8 is a system of m simultaneous equations in 



C i " * b i (s)a j (s)ds 

Aij = b ± (s)a (s)ds 



C i — A .2 A, .C . 
1 J=1 J-O J 



A2.6 

A2.7 

A2.8 



the m unknown constants C^. 

Equation A2.8 now represents a matrix eigenvalue problem. This 
can be seen by expanding A2.8 and writing it in the form of A2.9. 
where a = l/\. Equation A2.9 can then be written in the form of 



0- An 
” A 21 



0-1 

\ 



“*12 “ A lm 



J_A 22 



-A, 



2m 



-A. 



m2 



CJ-A, 



mm 





C 1 




“o“ 




C 2 


a 


0 




• 




• 




• 




• 




• 




• 




_ C m_ 




0 



A2.ll, which is a matrix eigenvalue problem, 
(al - A)C - 0 



A2.9 



A2.10 



A2.ll 
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The matrix eigenvalue problem A2.ll must now be solved. The 
first requirement for solutions to exist is that the determinant of 
the matrix be zero as in A2.12, This specifies an m*'- order 
polynomial in 0 which when solved gives the m - real roots or 
eigenvalues 0^j i - l,m. For each of these 0 ^ there must exist an 
eigenvector such that A2.13 is satisfied. 

det (01 -A) =0 A2.12 

(0 i I-A)E i . 0 A2.13 

Once these eigenvalues and eigenvectors are known, the first 
solution can be expressed by A2.14. It must be emphasized that there 
are m - solutions W^(z) to this problem, one for each eigenvalue. 

01 } W 1 (z) fn A2.14 

Thus, the eigenfunction solutions to the equation A2.1 axe given 
in A2.15. 

W^a) -1\ a (z)E A2.15 

1 j»l 1 J ij 
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APPENDIX A3 



SOLUTION OF A HOMOGENEOUS FREDHOLM II-INTEGRAL EQUATION BY CONVERSION 
TO LINEAR ALGEBRAIC EQUATIONS 

This appendix presents an approximate solution technique for 
determining the eigenvalues and eigenfunctions of an integral equation 
A3.1 by dividing the integration interval into equally spaced in- 
crements and changing the integral to an algebraic summation. The 

W(z) =» \ / C K(z,s)w(s)ds A3.1 

o 

integral equation eigenvalue problem is thus changed to a matrix 
eigenvalue problem. The approximate solution W(z) can then be 
expressed in terras of the eigenvalues and eigenvectors of the matrix 
eigenvalue problem. References for this material are: (H-13) (L-8) 

(P-4) (G-7) (S-ll) (T-3) (V-4) (S-14) (D-5) and (M-6). 

The first step in this solution technique is to divide the 
interval [0,C]| into n equal parts of length 6 as in equation A3. 2. 

6 = C A3. 2 

n 

Next, the functions W(z) and K(z,s) can be evaluated along the 
interval as A3. 5 and A3. 6 at each of the points given by A3. 3 and A3. 4. 



K(i6,j6) - 


(i,j = 1 , 2, •••, n) 


A3. 3 


W(i6) . W i 


(i « 1, 2, •••, n) 


A3. 4 


z ^ *» i6 


(i = 1, 2, ••*, n) 


A3. 3 


S, » j6 


(j = 1 , 2, n) 


A3. 6 



With these relationships established, the integral in A3.1 is 
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approximated by a summation over the interval £o,c] represented by 
the indices £l,n]] as in equation A3. 7. 

n 



W, 



X 2 
J-l 






*3.7 



Equation A3. 7 represents a matrix eigenvalue problem similar 
to equation A2.8 in Appendix A2. This equation A3. 7 can be placed 
in eigenvalue form A3. 9 by defining 0 as in A3. 8* similar to equations 
A2.9» A2.10, and A2.ll in Appendix A2, 

0 - 1_ A3. 8 

X 6 

(el-A)W =0 A3. 9 



Once the eigenvalues 0, and eigenvectors for equation A3. 9 
axe found, the approximate solution to equation A3.1 can be expressed 
as A3. 10. 

Wj.(z) = K(z,j 6 )E i j A3. 10 
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APPENDIX A4 



DIGITAL COMPUTER PROGRAM AND OUTPUT FOR EVALUATION OF PROBLEMS IN 
CHAPTER L4 

This appendix presents a listing of the Fortran IV statements 
in the program and subroutines used to calculate the numerical values 
of the four analytical solutions presented in Chapter L4, The numeri- 
cal results are first listed, then plotted, and finally printed out in 
matrix form. Tests made with this program showed that in all cases 
less than ten terms in the series were used for three-decimal-place 
accuracy. 

There are three subroutines used in the program: MXOUT, PLOT, 

and LOC. Subroutine LOC was taken, less comment cards, directly from 
reference (l-3) , the Scientific Subroutine Package, Subroutines MXOUT 
and PLOT, used here, represent significantly modified versions of the 
MXOUT and PLOT subroutines described in reference (1-3). The main 
program and subroutines are written in Fortran IV, and all runs were 
compiled and executed using the WATFOR compiler on the IBM 360/65 
computer at the M.I.T. Computation Center. The entire program 
compiles in 5»8 seconds and executes in 5«&5 seconds. 
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C PROGRAM TOT CALCULATE 4 TRANSIENT RE S PONSE S * UA l , UA2 1 UB 1 ,UC 1 

C PROGRAM TO CALCULATE 4 ANALYTIC SOLUTIONS AS A TEST 

DIMENSION UA( 100 J ,UAT ( 100 ) ,UB ( 100 ) , UC ( 100) 

__1_0_P_ fOR M A T (_3_X_ t J Xf 1»3 X_,J J= J x 7 X ,_*_U_A 0 i NE = !_t_6Xj_* U AT W 0= ' '_t_6_Xj _ 

1 * UBON E= • »6X»'UC0NE='»//> 

L°JL f OR M AT (_1_X j_* U A I W 0 _NOJ_ 1 _ P E RC_E_NJ_ A C CURAT E X = •_ t_F.7_._3j _ 

1' T=»,F7.3> 

103 F0RMAT(1X, »UAONE NOT ONE PERCENT ACC URATE , X= • , F7. 3 , 

_______ 

_ _l_C4_ F CR M A T (_1_X j J U B ON E U CONE^ _N_OJ_ ONE _PER_CENI _ A CCUR A T E x X= 

1F7.3,*T=* » F 7. 3 ) 

i_0_5_ FORMAT <_1_X_, 2 F 5 2 ,_3_X_ f 4 E 1 2 5 > 

PI=3. 14159265 

DELTT=0.025 

DELTX=0. 1 

MJN=Lt.OE-OJ 

A= 1 . 1 1 8 

Bf PI*PI_*_A*A 

BN=0 .25*8 

WRITEI6, 130 ) 

C THIS SECTION IS TO SET UP OUTPUT PLOTS AT 0.1 INTERVALS 

C__IN_X. 

XP=0.0 

D0_20_0_ 1=1,10 

X P=X P+DELTX 

UA I I ) =XP 

UAT ( I )=XP 

UBillj^XP 

UC ( I ) = XP 

__2J^_CCNTI_NU_E 

1=11 

T= C , 0 

DO 10 J=1 »9 

T=I+DE_L_TT 

X=0 .0 

DC_20__K=_l_,i0 

X=X+DELTX 

C TIME-SPACE GRID IS ,\’CU SET 

C STEADY STATE CALCULATIONS 

UAGST=_X 

UATST=1 .0 

AC = K_P./_( _A*A1 

UBOST = ( 1.0-EXP( AC*X> )/ll .0 -EXP (AC) ) 

UC CST=(1. O-EXP (- AC*X ) )/( 1 . 0-E X P ( - AC ) ) 

C TRANSIENT CALCULATIONS 
SN=-1 .0 
UACTR=0. 0 

_ UATTR = 0.0 _ 

UBCTR=C . 0 

ATB=2.0*EXP(( — 0 .25*T+0.5*X-0.5 )*AC) 

ATC=2.0*EXP( (-0.25*T+0. 5-0 . 5*X ) *AC ) 

JAOfO 

I AT=0 
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I B0=0 

IFJN=0 

C USE CMY THE FIRST 400 TERMS OR LESS OF THE SERIES 
DO 30 N=l,4C0 
IFIIFIN- 3)53, 54,54 

53 RN-FLCAT(N) 

AN=RN*P I 

If i I A ' T I_4_I_, 41 _, 4 2 _ 

41 RTN-FLOAT (2*N-1 ) 

ATN=RTN*PI 

CCSN = COSI ATN*X*0 .5 ) 

e _xpbn=_e_x p_i -r IN ?.R ' tmj *b n i _ 

UATT=UATTR 

U A T TR= IMTJR + (S. N *COS N * E X P B N )_/J AJN *0 . 5_)_ 

IFIABSI UATTR)-MIN)44,44,6l 
61 IF (ABS<(UATTR-UATT)/UATTR)-0.0 1 >43,43*42 

44 UA TTR=0. 0 
_ __4_3 __IAT = 1 

IFIN=IFIN+1 
_42_ N= S I_NJ AN * X I 

EXPB= EXP I— RN*RN*T *B ) 

I F ( I A0)45*45* 46 

45 U AO=U AOTR 

_UAG IR=JJ_AOJ R + IS N > *. S I N ME X P B )_/_AN 

IFIABSI UACTR) -MIN >48,48,6 2 
_ 62_ IFUBSUUA 0 T R - U A0_)_/ U AOTR L-0_._01 ) 4 7 ,47,4 6 

48 UA0TR=0. 0 

47 I A0= 1 

IF IN= IFIN4-1 

_4 6_ I F i IB 0 )_49_, 4 9 L 5C 

49 U80-=UB0TR 

SUM S Q = 0_._2 5 * A C* AC+A N *_A N 

UBOTR=UBOTR+I SN*AN*S INMEXPB ) /SUMSQ 

IF(ABSIUBQTR)-MIK352,52?6 3 

63 IF (ABSI IUBOTR-UBO)/UBOTR >-0.01 >51,51,50 

52_yBOTR_^0_.0 

51 I E0=1 

ifIN=LF_I_N+l 

50 SN=-SN 

30 CONTINUE 

IFIIAT)55,55,56 

_5 5_ W RII E_I_6_,i 02 _)_X t_J_ 

56 I F ( I AO ) 57 , 57, 5 8 

_57_ W R IT E_I_6_, J 0 3 J_X $_T_ 

58 IF{ 180)59,59,54 

59 WR I TE ( 6 , 1 04 ) X , T 

54 ATR-2 . 0*UAOTR 

U_A ON E = U . ADST + ATR 

ATTR=2 ,0*UATTR 

U A J W 0= UATST + AJIB. 

BTR=ATB*UBOTR 

CTR=ATC*UBOTR 

UBONE=UBOST+BTR 

-JJCONE =UCOSJ+CTR 

C SAVE VALUES FOR PLOTTING 
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UM I ) =UAONE 

U A Tj. I )_=UATWO_ 

UB ( I ) =U80NE 

UC ( I ) =UCONE 

1 = 1*1 

WRI T E L6 »_1_0 5_) X L TfJUAONE _, U ATwQ_,JJ 8CNE L UCO_NE_ 

20 CONTINUE 

1_C_ CONTINUE- 

C PLOT COMMANDS 

C MATRIX PRINTOUT 

N0= 1 

MP = 10 

NP = 10 

ms=o 

L I N S= A 5 

I PQS=65 

ISP=1 

_NL=46 

NS = 1 

CALL _PL_CJ_( NO _U A ,_N_P_, PP ^NL »_N_S_) 

CALL MXOUTINO, U A ,NP , MP , MS , L INS , I POS, I SP > 

NONO+1 

CALL PLCTINO,UAT,NP,MP,NL,NS) 

CALL _MX_0_UTj{ N 0 »_U A T_,_NP, MP t M S ,_L_I N I S , IPO: S_,_I S P J 

NC=N0+ 1 

CALL _P LCJJ NO , _U B ,_NP_, M P ± N L ,_N_SJ 

CALL MXOUTINO, UB ,NP , MP , MS, L IN S, I POS , I SP ) 

N0=NQ*1 

CALL PLCTINO, UC , NP , MP ,NL , NS ) 

CALL _M_X_GUJi NO i. _UC_,_NP_, M P t MS_,_L_INS_, I P 0 S_,_I S P 1 

CALL EXIT 

END 
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SlB'RCfOT I N'E PLOT ( NO , A , N , M , NL , NS ) 

C ****±±^*J?***±±* ^I_CH_AEL_N .HAYES _ * * [** * * * * 

DIMENSION OUT ( 10 I ) , Y PR( 1 1 ) , I ANG ( 9 } , A { I> 

INTEGER IDUM/ * 1 «/, I ANG/ » 1 ♦ , » 2* , * 3 * , * 4 * , « 5 » , * 6 * , * 7 » , 

1 *8 * » *9* / 

JNJEGER _0_UJ 

C THESE LIMITED FORMATS ARE FOR 60 SPACE PRINTOUTS FOR 

C J H E S I_S JJ_S_E_. 

1 FORMAT( 1H1,28X,7H CHART ,13,//) 

2 FORM AT ( 1H , F8 . 3 , IX ,« *■ , 51 A 1 , * * « ) 

3 FORMAT ( IH ) 

7_f ORM AT (_I_H__, J.OX » 36H_*_ _ * * * * * * 

I15H * * *) 

_8_ f 0 RM AT (_1_H0_, 8 X , 1 IF 5. 2) 

C NTH IS NUMBER CF SPACES TO BE USED, EITHER 101 OR 51 

N TH=5 1 

NLL=NL 

I f i N S )_1_6,1 6 1 0 

10 DO 15 1=1, N 

J50_14_Jj=_I_,N 

IF(A(J)-A(I 1)14,14,11 

11 L= I-N 

LL= J-N 

00_12 _K=_1_,M 

L=L + N 

LL=LLt/l 

F= A ( L ) 

A(L) = A(LU 

12 A(LU=F 

14_C0NTI_IW_E 

15 CONTINUE 

1 6_ I F iNLLJ 20_,i 8i2C 

18 NLL= 50 

20 WRITE(6«1>N0 

WRITE(6,7) 

BLANK=_0 

XSCAL=( A(N)-A( 1 ) )/( FLOAT (NLL-1 ) ) 

M!=N+1 

YM AX=— 1 .0E75 

YM I N= 1 » 0E75 

M2=M*N 

D0_40_ _J= M_i_, M2 

IF ( A ( J ) .GT. YMAX > YMAX=A(J) 

I F _ (_A <_J)_ ^ L T . Y_M_I_NJ_ YM LN= MJJ 

40 CONTINUE ‘ 

YSCAL= CYMAX-YXIX 3/5C.C 

X6= A ( 1 ) 

L=_l i 

M Y=M— 1 i 



45 F= I -1 

XPR=XB+F»XSCAL 

IF (XPR-A(L) ) 50 ,50 ,70 

5 C_ D 0_ 5 5 _ I_Xf_l_,_NJH 

55 OUT ( IX ) =BLANK 
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DO 60 J= 1 » MY 

JAf L±J*N_ 

JP=(( A(LL)-YMIN)/YSCAL)*1.0 

OUT (JP)=IANG(J) 

60 CONTINUE 

C .PRO PER FORM AT _ C A R_C_S_ MUST BE. JJS E D 

WRITE(6,2)XPR, (OUT (I 2) »IZ=1,NTH) 

L=L+ 1 

GO TO 80 

70 WRITE(6«3> 

8C 1 = 14-1 

If I I-NLU45* 84 L e6_ 

84 X PR=A ( N) 

GO_TO _50 

86 WRITE(6,7> 

YPR ( 1 ) =YMI N 

DO 90 KN= 1 » 9 

___9_0_ Y P R 1KN+1J = Y P R JK_ N )_+Y S C A l * 5 ,_0 

YPR( 1 1 ) =YMAX 

WR_ITE L6_,_8J_( Y P R (_ LP_b_I £= I ± LU 

RETURN 

END 



- 194 - 



SUBROUTINE MXOUT (I CODE , A , N , M , MS , L I NS , I POS, I SP » 

C Jf*****^*^**** Jili-HAEL N j__HAY_ES ^*J?J?***fl*^*J?J?**** 

DIMENSION A ( 1 > t B(8) 

I FORMAT! 1HO, 5X» 7HMATRIX ,I5,6X,I3,5H ROWS, 6X, 13, 

18H COLUMNS,/) 

2_ FOR MAT(_1_2_X_, 8HC OL UMN_ 3X_,_I_3_,Jl 0 X ) ) 

3 FORMAT ( 1H ) 

A_ FOR MAT _, 7 X 1 4 H ROJ»_ _I 3 (_EJ_6_.6ii 

5 FO RMATII HO ,7X, ARROW , 13, 7IE16.6) ) 

6 FCRMAT { 1 H 1 > 

7 FORMAT (16X,13HST0RAGE MODE , 1 1 , 7X , 7HGROUP ,12,/) 

WRJ[TE (_6_,_6) 

J=1 

N EN D= I_P_0_S_/_1 6- 1 

LEND=(LINS/ISP)-2 

I PAGE = I 

10 L S TRT = 1 

_20_ W R 1 1 E (_6_,_1_) JCODE ,_N_,_M 

WRITE(6,7)MS,IFAGE 

JNT=J + NE ND-JL 

I PAGE= IPAGE+1 

31 I F I JN T-M ) 33 , 33 ,3 2 

32 JNT=M 

33C0NTLNUE 

WRITE(6,2 )( JCUR,JCUR=J,JNT) 

SP-JJ _35_, 3 5 j_Ai p 

35 WR IT E( 6, 3 1 

AC L TEND^L STRT+L END-1 

DO 80 L=LSTRT»LTEND 

DO 5 5 _K=_1_,NEND 

KK=K 

JIzJtK-J 

CALL LCC(L,JT»IJNT ,N , M , MS ) 

B( K?=0 .0 

I F I IJNT)5C,50,A5 

_A_5_BI K i= A (J JNJi 

50 CONTINUE 

I F i JT ■ ^MJ 5 5jt 6 C ,_60_ 

55 CONTINUE 

60 IF(ISP-l) 65,65,70 

65 WRITE(6,A)L,(B(JW),JW=1,KK) 

GO_TO _7_5 

70 WRITE(6,5)L,(B( JW), JW=1,KK) 

_75_ I F IN- LJ_8 5 _,85 L 8_C_ 

8C CONTINUE 

LSTRT-LSTRT»LE.\D 

GO TO 20 

_8_5_ _IF_( JT- M )_90_,95_,95 

90 J = JT + 1 

gc_jo_lq_ 

95 RETURN 

END 
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-SJJ B6.CLIiTJLN£_ L QC L U J.» J B * N 

C jot***#****##**#*#*#******##**#**##*#**** **#$***##*** 

I X = I 

JX = J 

IE LMS-_lJ_ IQ j. 2.0 13_0 

10 IRX=N*( JX-1)+I X 

GQ_U1_3_6 

20 IF(IX-JX) 22,24,24 

22 IRX=IX+(JX*JX-JX )/? 

GO TO 36 

2j4_IR X= JJL+J J 1 * I X.- 1 XJ7 2 

GO TO 36 

3JD_IRX=<1 

IF(IX-JX) 36,32,36 

32 1 RX= IX 





36 I R = I R X 
RETURN 
END 
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x= 


T= 


UAONE= 




U AT W0= 




U50N£= 








C.10 


0.02 


0.3059 5E- 


-03 


0 . 37324E' 


-03 


0.21213E— 03 


0 . 43643E • 


-C 3 


0.20 


0.02 


0. 13734E- 


-02 


0 .3569 IE- 


-02 


0.993C1E- 


-03 


0. 18837E- 


-02 


0.30 


0.02 


0.50529E- 


-02 


0. 49 80 IE- 


-02 


0.38026E- 


-02 


0 . 66583 E- 


-0 2’ 


0.40 


0.02 


0. 16 39 IE- 


-01 


0.166C5E 


-01 


0. 12842E- 


-01 


C.20754E- 


-01 


0.5C 


0.02 


0. 32329E- 


-01 


0. 45695E- 


-01 


0 . 26402E- 


-01 


0 . 39390E- 


-0 1 


0.60 


0.02 


0. 10959E 


00 


0. 10548E 


00 


0 .93069E- 


-01 


0.12817E 


00 


0.70 


0.02 


0.22939E 


00 


0.22987E 


OC 


0.20286E 


00 


0. 25789E 


00 


0.80 


0 .02 


0.4237CE 


00 


C.42378E 


00 


0 .39023E 


00 


0.45795E 


00 


0.90 


0.02 


0 . 689 2 IE 


00 


0.68916E 


00 


0.66130E 


00 


0. 71638E 


CO 


1.00 


0.02 


0. 1C000E 


01 


C. 10000E 


01 


C. 10000E 


01 


0. 10000E 


Cl 


0.10 


0.05 


0. 9037 IE- 


-02 


C.12726E- 


-01 


0. 62542 E- 
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MATRI X 


3 


10 


ROWS 


10 COLUMNS 










STORAGE MODE 


0 


GROUP , 


1 








COLUMN 


1 




2 




3 




ROW 


I 




0. 100000E 


01 


0. 1 00000 E 


0 1 


0 . 1000CCE 


01 


ROW 


2 




0. 900000 E 


00 


0. 661296E 


00 


0.745336E 


CO 


ROW 


3 




0.800000E 


oc 


0. 390233E 


00 


0. 525795E 


00 


ROW 


4 




0 • 700000E 


00 


0. 202861 E 


00 


0. 349673E 


CO 


ROW 


5 




0. 600000E 


00 


0 . 930 687E- 


-01 


0.218502E 


CO 


ROW 


6 




0 . 500000 E 


00 


0.264024E- 


-01 


0. 127269E 


CO 


ROW 


7 




0 • 400000E 


oc 


0. 128419E- 


-0 1 


0.699988E- 


-Cl 


ROW 


3 




0. 300000E 


00 


0 . 380260E- 


-02 


0. 356159E- 


-01 


ROW 


9 




0.2 COOCOE 


00 


0.993013E- 


-03 


0. 165389E- 


-01 


ROW 


10 




0. 100000E 


00 


0.21213 3E- 


-03 


0.6 2541 6E- 


-C2 


MATRI X 




3 


10 1 


ROWS 


10 COLUMNS 










STORAGE MODE 


0 


GROUP , 


2 








COLUMN 


4 




5 




6 




ROW 


1 




0 . 1 COOOOE 


01 


0. 100000E 


0 1 


0. lOOCCOE 


01 


ROW 


2 




0.783344E 


oc 


0.806019E 


00 


0.821203E 


00 


ROW 


3 




0. 591927E 


00 


0.632494E 


00 


0. 659998E 


CO 


ROW 


4 




0 . 430447E 


00 


0.482313E 


00 


0. 5181 77E 


00 


ROW 


5 




0. 300572E 


00 


0. 356637E 


oc 


0. 396426E 


00 


ROW 


6 




0 . 200998E 


cc 


0. 254890E 


00 


C.294283E 


00 


ROW 


7 




0. 128274E 


00 


0.174963E 


oc 


0.210224E 


00 


ROW 


8 




0. 772659E- 


-01 


0. 113586E 


cc 


0. 141S65E 


00 


ROW 


9 




0. 424448E- 


-01 


0 . 667886 E- 


-01 


0.862 302E- 


-Cl 


ROW 


10 




C. 183465E- 


-01 


0.3031 76E- 


-01 


0.400413E- 


-Cl 


MATRI X 




-a 


10 ROWS 


10 COLUMNS 










STORAGE MODE 


0 


GROUP , 


3 






COLUMN 


7 




8 




9 




ROW 


1 




0. 100000E 


01 


0. 100000E 


0 1 


0. 100000E 


Cl 


ROW 


2 




0.831905E 


oc 


0.839613E 


00 


0.845214E 


00 


ROW 


3 




0. 6 7948 5 E 


oc 


0.693550E 


00 


0.70378CE 


CO 


ROW 


4 




0.543 800E 


00 


0 . 562358E 


oc 


0.575874E 


00 


ROW 


5 




0,4251 59 E 


00 


0. 446060E 


00 


0.461308E 


00 


ROW 


6 




0.323077E 


00 


0.344124E 


00 


0 . 3595CSE 


CO 


ROW 


7 




0. 236322E 


00 


0. 2 5 5493 E 


oc 


0.269533E 


00 


ROW 


8 




0 • 1 63038E 


oc 


0. 178660E 


00 


C. 190122E 


00 


ROW 


9 




0. 100923E 


00 


0. 1 11804E 


00 


0. 119797E 


00 


ROW 


10 




0. 474349E- 


-01 


0. 529225E- 


■Cl 


0 . 569578E- 


•01 


MATRI X 




3 


10 ROWS 


10 COLUMNS 
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STORAGE MODE 0 



GROUP , 4 



{ 



COLUMN 10 



ROW 1 C.1COOOOE 01 

ROW 2 0 j.84929 8E_ C C 

ROW 3 0. 711242 E 00 

ROW 4 0 .5 8_5_7 3 9 E _ 00 

ROW 5 0. 472445E 00 

ROW 6 C «370754E CO 

ROW 7 0. 279803E 00 

ROW _8 0 •_ 1_ < 58_5J. 2 E _ 0 0 

ROW 9 0. i 25652E 00 

ROW 10 C . 599143E-0 1 
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— 


* 


* 


* 
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MATRI X 



4 



10 ROWS 



10 COLUMNS 







STORAGE MODE 


0 


GROUP , 


1 








COLUMN 


1 




2 




3 




ROW 


l 




0.100000E 


01 


0. 100000E 


01 


C. 100000E 


01 


ROW 


2 




0 . 900000E 


oc 


0. 716376E 


00 


0.807416E 


00 


ROW 


3 




0. 800000E 


00 


0.457947E 


OC 


0.617032E 


00 


ROW 


4 




0 . 700000E 


oc 


0.257891E 


00 


0.444528E 


00 


ROW 


5 




0.600000 E 


00 


0. 128 171 E 


00 


0.30091 IE 


00 


ROW 


6 




0 * 500000E 


cc 


0.393896E- 


-01 


0. 1 89869E 


00 


RCW 


7 




0. 400000 E 


00 


0.207545E- 


-01 


0. 113127E 


CO 


ROW 


8 




0. 300000E 


00 


0.665832E- 


- C 2 


0.623550E- 


-01 


ROW 


9 




0. 200000E 


oc 


0. 188375E- 


-02 


0.3 1 367 3E- 


-01 


ROW 


10 




C. 100000E 


00 


0 . 436425E- 


-03 


0. 128500E- 


-01 


MATRI X 




4 


10 i 


rows 


10 COLUMNS 










STORAGE MODE 


0 


GROUP , 


2 








COLUMN 


4 




5 




6 




ROW 


1 




0 • 100000E 


01 


0. 100000E 


0 1 


C.IOCOOOE 


01 


ROW 


2 




0. 848590E 


oc 


0.8731 54E 


00 


0 . 88960 3E 


CO 


ROW 


3 




0.694639E 


00 


0.742246E 


00 


0. 774522E 


00 


ROW 


4 




0. 547213E 


00 


C.613148E 


00 


0.658741E 


00 


ROW 


5 




0.4 13934E 


00 


0.491 145E 


00 


0. 545941E 


00 


ROW 


6 




0. 299862E 


00 


0.380262E 


00 


0.439030E 


00 


RCW 


7 




0.207307E 


oc 


0.282761E 


OC 


0. 339748E 


CO 


ROW 


8 




0. 135273E 


oc 


0. 198860E 


00 


0 . 248368E 


00 


ROW 


9 




0 . 804995E- 


■01 


0.1 26669E 


oc 


0. 163541E 


00 


ROW 


10 




0 . 376939E- 


-01 


0. 622888 E- 


•01 


0 * 822664E- 


-0 1 


MATRI X 




4 


10 ROWS 


10 COLUMNS 










STORAGE MODE 


0 


GROUP , 


3 






COLUMN 


7 




8 




9 




ROW 


1 




0. 100000E 


01 


0. 100000E 


0 1 


0. 100000E 


01 


ROW 


2 




0. 901 196 E 


00 


0 * 909546E 


00 


0.9156 13E 


00 


ROW 


3 




0. 797390E 


00 


0.813896E 


00 


0 . 8 2590 IE 


CO 


ROW 


4 




0. 6913 15 E 


00 


0. 7149C8E 


00 


0. 732090E 


00 


ROW 


5 




0.585510E 


00 


0.6 14294E 


00 


0. 635293E 


00 


ROW 


6 




0 .481987c 


00 


0. 5 1 3387E 


00 


0. 536338E 


00 


ROW 


7 




0 • 381926E 


CO 


0.4 12908E 


00 


C.435598E 


CO 


ROW 


8 




0 . 285437E 


oc 


0.312787E 


00 


0. 332854E 


CO 


ROW 


9 




0. 191407E 


00 


0 . 2 12042E 


oc 


0. 227202E 


00 


ROW 


10 




0.974566E- 


01 


0. 108731E 


00 


0. 117022E 


CO 



MATRIX 4 



10 COLUMNS 



10_ ROWS 
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STORAGE MODE 0 



GROUP 



t 



4 



COLUMN 10 

ROW 1 C . 1 OOOOOE 01 

ROW 2 0._92_00_3 8 E OC 

ROW 3 0. 834659E 00 

ROW J* 0 ^7 44_6_3_1 E _ 00 

ROW 5 0. 65063 1 E OC 

ROW 6 0.553115E PC 

RCW 7 0. 452196 E 00 

ROW _8 0 .3475 4 3 E _gg 

ROW 9 0.238306E OC 

ROW 10 0 . 123096E 00 



COMPILE T I ME= 5.80 S EC , E XECUT I ON T I ME = 5.65 SEC 
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APPENDIX A5 



SUBROUTINES USEFUL FOR COMPUTING LPCM STEADY-STATE SOLUTION 

For completeness and for the convenience of the reader, this 
appendix present's a listing of the Fortran IV statements in the sub- 
routines useful for computing the LPCM steady-state response. As 
mentioned previously it is expected that there are easier and more 
efficient methods which could be developed and used to find the 
steady-state for special cases of the LPCM. These subroutines will 
solve the general LPCM steady-state with the subroutine AFCT given 
here. Several of the subroutines here are also used directly in PLPBV, 
namely EIGEN and GELG. 

Four of the subroutines used here are taken, less comment cards, 
directly from reference (l-3) * the Scientific Subroutine Package. 

These are subroutines LBVP, EIGEN, GELG, and LOC. The remaining 
three subroutines, AFCT, FCT, DFCT, were written specifically for 
the LPCM. One additional subroutine, OUTP, is required for operation 
of the subroutine LBVP and must be furnished by the user. All of 
these subroutines compile under WATFOR in a total time under 20 seconds. 
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0 SUB ROUTINE LBVP(PRMT,R,C,R, Y, DERY,NDIM, T HL F , A FCT , FCT , DFC T , O'JTP , 
1AUX,A) 

D I MENS I ON PRMT (1 ) , B( 1 ) ,C (1 ) , R ( 1 1 ,Y (1 ) , DER Y ( 1 ) , AUX ( 29 , 1 ) , A ( 1 ) 
IF(PRMT(3)*(PRMT(?)-PRMT( 111)2,1,3 

1 I HI F= 1 ? 

RETURN 

2 IHLF=13 
RETURN 

3 KK=-ND T M 
I B=0 

IC = 0 

DO 7 K = 1 1 ND TM 
Al)X(l S,K)=DERY(K) 

AUX( 1 ,K)=1 . 

AUX ( l 7 »K ) = 1 • 

KK =KK+ND I M 
DO 4 1 = 1 f ND TM 

T T = KK + 1 

T F ( B ( TI))S t 4, 5 

4 CONTINUE 

T R = I R + 1 

AUX ( l , K 1 =0 . 

5 DO 6 T = 1 , ND I M 
I T =KK+ I 

T F ( C ( m 17,6,7 

6 CONTINUE 
IC=IC+l 

AUX ( 17,K)=0. 

7 CHNTINUF 

IF{ TC-TR18, 11,11 

B H= PRMT { 2 ) 

PRMT ( 7 )=PRMT ( l ) 

PRM T ( 1 T = H 
PR MT ( R ) =-PRMT ( 3 ) 
on 1 I =1 1 ND I M 
9 AUX ( 1 7 , I) = AUX ( 1 , I ) 

I I =ND I M*NO I M 
DO in 1=1,11 
H=R( I 1 
B(T)=C(n 
in C( T)=H 

11 X= PRM T (2 ) 

CALL FCT(X,Y) 

CALL DFCTIX ,DERY ) 

DO 12 1 = 1 » NO I M 
AUX (18,! ) = Y (T 1 

12 AUXI 19,T)=DERY( T 1 
K= 0 

KK=C 

100 K = K + 1 

I F ( AUXI 1 7 , K 1 )10R, 108,101 

101 X= PRMT ( 2 ) 

CALL AFCT(X.A) 

SUM=0. 

GL=AUX( 18, K) 
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DGl=AUX( 19, K) 

T I =K 

DO 104 I = 1 » ND I M 
H=- A{ T I ) 

DERYU )=H 

AI.JX ( 2 r , I ) =R ( I ) 

Y{ I )=o. 

IF( 1-K)101,10Z,1 03 
io? v( n = i . 

103 DGL = OGL+H-*AUX( IP , I) 

104 I T = I I + N0 1 M 
XF ND = PRM T ( 1 1 
H=.0625*( XEND-X) 

I sw=n 

GOTO AO^ 

P5 I«=(THLF-1^> 106,106,117 

106 DO 107 1=1 , NDIM 
KK=KK+ 1 

H=C( KK ) 

R ( I)=AUX ( ?C , I ) +H*SUM 
I I =1 

DO 107 J=1 , ND T M 

B( rn = B< m +h*y( j > 

107 11=11+ NDIM 
GOTO 1 09 

108 KK = KK + ND I 

loq IF(K-N0IM)100,110,U0 
no X= PRM T ( 4 ) 

CALL GELG(R,B,NDIW,1 ,X, I ) 
IFCDlll ,112,11? 

111 I HLF=1 4 
RETURN 

112 PR MT ( 8 ) =o . 

I HLF=- I 
X=PRMT (1 ) 

XFND=PRMT{ 2> 

H=PRMT(?) 

DO 113 1=1, NDIM 

113 Y(I)=Rm 
ISW = 1 

114 I SW2= 1 2 
GOTO ?00 

116 I SW3 = - l 
GOTO 300 

116 1F( IHl F)4O0, 400,11 7 

117 RF TURN 

200 CALL AFCT(X,A1 

IF(IRW)?oi,?01,?O e > 

?01 LL=0 

nn ?03 m=i,ndim 
HS = 0. 

DO 207 L = 1 , NOI M 

LL=LL+ 1 

?C 7 HS=HS-A(LL)*Y(L) 

203 OERY{M)=HS 
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GOTO ( 802, 804, 806,407, 418, 41 8,608 ,617 ,632, 634,421 ,118),! SW 1 

205 CALL FCT(X,DERY) 

DO 207 M = 1 , NO I M 
LL=M-NDIM 

H$ =0 • 

on 206 L = 1 , NO 1 M 
LL=LL+NOTM 

206 HS=HS + 4 ( LL ) *Y{ L) 

207 OERY{M)=HS+OERY(M) 
onm 204 

800 I«= ( I SW)301 ,301 ,308 
701 CALL FCT(X,R) 

GU«0. 

0GU=0 . 

on 302 L = 1 , NOIM 
GU=GU+Y(L) *R(L) 

30? DGU=DGU+DERY( L >*R(L ) 

CALL OFC T ( X , R ) 

On 3no L=1 , NO I M 
3 07 DGU= OG U+Y ( L ) *R ( 1 ) 

SUM=SUM+.5*H*( (GL+GU)-*. 1666667*H*( OGL-OGU) ) 

GL = GIJ 
DGL =DGU 

304 !<=( ISW3) 116,422, 61 8 

7H3 CALL nUTP(X,Y,DERY,lHLF,NOIM,PRMT) 

IF{PRMT( 8) ) 11 7,304, 117 

400 N= 1 
XGT=X 
IHLF=p 

Dn 401 1 = 1 , NOT Y 

AUX ( 16, T )=0 . 

AUX ( 1 ,I)=Y{1 ) 

401 AUXT8 , I)=OERY{ I) 

I GW1 = 1 

GOTO 8 on 

402 X=X+H 

DO 4"7 I-l,NniM 

403 AUX ( 2 , I ) =Y ( I ) 

4° 4 I HLF= THLF+1 

X =X-H 

On 408 1=1, NOIM 
4ns AUX{4, T)=AUX( 2, T ) 

H=.6*H 
N= 1 

I8W1=7 
GOTO 80° 

406 X=X+H 

I SW2 = 4 
GOTO 200 

407 N= 2 

DO 408 1=1, NOIM 
AM X ( 2, T)=Y( I ) 

408 AUX (0, T ) = 0E° Y ( I ) 

T 5W1 =8 

GOTH 8 00 
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409 Da 414 1 = 1 , NO I M 

7 = AR $ ( Y ( I) ) 

1 F ( Z - 1 . >410,411 ,411 

410 Z=l. 

411 DFLT= . 06666667* A8G(Y( I )- AUX (4,1)) 

IF (T SW) 417,413 ,412 

412 DFl T=AIJX(15,I )*DEl T 

413 I F (DFl T-Z*PRMT (4 ) ) 41 4, 41 4, 420 

414 CONTINUE 
X = X + H 

I SW2 = 5 
GOTO 200 

41 5 DO 414 1=1, ND I M 
AIJX (3 , T )=Y (II 

416 AUXI 10,1 ) =OERY( I ) 

N=3 

I SW1 =4 
GOTO 5 nr> 

417 N= 1 
X= X+H 

I SW2= 4 
GOTO 200 

41 8 X=X$T 

DO 410 1=1, ND T m 
AUXI 11 , T 1=DEPY( I ) 

4 1 00 Y ( T ) = AIJX (1 , I )+H*( .37 5* AUX (8,1 )♦. 70 16667* AUX ( 0,1) 
1-. 2083333* AUX I 10, I) + .0 41 66667*DERY ( I ) ) 

420 X=X+H 
N=N+1 
TSW2 = 1 1 
GOTO 20 0 

421 TSW3=o 
GOTO 300 

422 TF (N-4) 423, 600,600 

423 DO 424 T=1 , ND I M 
AUXIN, n=V( T ) 

404 AUXI N+7 t I) =DE RY I I ) 

IF (N — 3 14 25, 427, 60 0 
425 DO 424 1=1 , ND I M 

DEL T= MJX I 0 , I ) +AUXI 0,1 ) 

DEt.T=DEl T+DEI.T 

4 26 Y I I ) = AUX( 1 , T ) + .7333333*H*( AUX I 8 , I M-DEIT+ AUXI 10,1)) 
GOTH 42^ 

427 DA 428 I =1 , ND I M 

OFl T=AUX(O t n+AUX( 10,1 ) 

DFLT = DEl T+DFLT + DELT 

428 Y( I ) = AUXI 1 , I ) + .3 75*H*( AUXI 8, I ) +DFL T + AUX (1 1 , I) 1 
GOTO 42 0 

420 IF! IH1 F-10) 404 ,430,43A 

470 I HLF = 1 1 
X = X-*-H 

I F ( T S W ) 105, 105,114 
500 Z=X 

DO 501 1=1, NDI M 

X=H*AUX(N+7, 1) 
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AUX( 5, n = x 

6 A 1 Y ( I ) * AUX ( M t I ) + .4 *X 
X=Z+.4*H 
I 8W2=1 
GOTO 2 00 

5 A 2 DO 503 1 = 1 , NO I M 

X=H*OFRYU ) 

AUX( 6 ,!)=X 

503 Y( n = AUX(N,I) + .2969776*AUX(5,I ) +.1 58 7 595* X 
X=Z+.4557372*H 

I SW2 = 2 
GOTO 200 

504 Dn 505 I=1,NDIM 
X=H*OERY{ I ) 

AUX (7,1 ) =X 

50 c; Y ( I ) =AUX(N, I ) + . 2181 n 04*A!)X (6,1 )-3 . 0 50058* AUX ( 8 , I)+3. 832866* X 
X = 7+H 
I SW2=^ 

GOTO 200 

5 06 on 507 I=l t NOIM 

5 07 *y( i ) = AUX ( M, I ) +. 174760 3* AUX ( 5, I 1-.5514 807* AUX ( 6, I ) 

1 + 1 . 205536*AUX( 7,T)+.1711R48*H*0FRY( H 
X=Z 

GnT 0 ( 402 , 406 , 400,417 ), ISW 1 

600 I $TFP=3 

601 I F ( N-8 >604,607,604 

602 OH 603 N = 7,7 

Dn 600 1 = 1* NO I M 
AUX ( N-l , I )=AUX(N» I ) 

603 AUX ( N + 6 , I ) =AU X { N+7 , T ) 

N=7 

60 4 N= N+ 1 

DO 60 6 I = 1 , NO I M 
AUX(N-l,T)=v( n 

605 AUXlN+6, I)=OERY( I ) 

X= X+H 

606 ISTFP= ISTEP+1 
On 607 I=1,N0TM 

ODFLT = AUX( N-4, T } + 1 . 3 333 33 *H* ( AUX (N + 6 , I ) +AUX(N + 6 , I )-AUX(N+5, I ) + 
1 AUX ( N + 4, I )+AUX( N+4, I ) ) 

Y ( T)=DELT-. 925619 8* AUX (16,1) 

607 AUX ( 1 6, T ) =DFLT 
I SW2 =7 

GOTO 200 

608 00 600 I=1,NDIM 

0 0F1 T= .1 2 5*( 9. * AUX (N-l t I ) - AUX (N - 3 , I ) +3 .*H*( OER Y ( I ) +AUX(N+ 6 , I 1 + 
1 AtJX(N+6,I>-AUX(N+5, I) ) ) 

AUX( 16, I)= AUX (16,1 1 -DE 1 .T 

609 Y { I ) =OELT+. 07438037* AUX(16,T ) 

OELT=o. 

no 616 1 = 1 , NO 1 M 
7= ARS ( Y( I ) ) 

IF (7-1. >610, 6 l 1,6 11 

610 7=1. 

611 Z=ARS(AUX(16,I 11/Z 
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IF(TSW)617, 613,612 
61? Z= AUX I 1 5 » I ) *Z 

613 TFI7-PRMTI4) >614,614,628 

614 TFIOEL T-Z >615, 616,616 

615 DFL T =Z 

616 CONTINUE 
I SW2 =8 
GOTO 200 

61? TSW3=1 

GOTO ?0O 

618 TF(H*(X-XENO> >619,621,671 
6 IP IF ( ABS ( X-XEND>— . I * AB S I H) ) 6 71,620,620 
670 IF ( DEL T-.n?*PRMT{ 4) ) 672,677, 601 
6?1 IF(I SW) 10S, 105,1 17 
62? I F ( IHLF ) 601 ,601 ,623 

627 TF(N-7)60l,624,6?4 

674 IF( TSTEP-4)60l,626,6?5 

675 IMOD=ISTED/? 

I F ( I STEP- 1 MOO- 1 MOD) 60 1 ,6 76,60] 

626 H= H+H 

thlf=ihlf-i 
T STFP =0 

00 627 T=1,NDIM 
AUXIN-1, T)=AUX(N-?, I) 

AUXIN- 2, I ) =AUX( N-4, I ) 

AUXIN-?, I )=AUX(M-6, I ) 

AUX(N+6,n=4UX(N + 6,l ) 

AUXIN + 5, 1 ) = AUX ( N+7, T ) 

AUXIN+4, I) = AUX( N + l, I ) 

OE I.T= AUX ( N+6 , I ) + AUX I N + 5 , I ) 

DFLT=DFLT+DELT+OELT 

62 70 AUX 1 16, 1 1=8. 0 629 6**1 VI I) -AUXIN-3, I) ) -3 . ?6 1 1 1 1 *H* I DFR Yt T ) +OEL T 
1 +AUXI N+4 , T ) ) 

GOTO 601 

628 I HL F= IHLF + 1 

IF! IHLF-10)630,67O,679 
679 IF! ISW)105,105,114 
6?0 H=.5*w 
TSTE P=0 

On 631 I=1,NDIN 

OYI I ) = . 0039Q 625*180. *AlJX( N-l , T ) + 1 7 5 . *AU X I N-? , T ) +40 . * AUX (N-3 , I ) ♦ 

1 AUX ( N-4, I))-. 11 71 8 75* ( AUX I N+6, T ) -6 .*AUX{ N+5, I)- AUX I N+4, I ) ) *H 
OAUXI N-4, T > = .00390625*1 1?.* AUX ( N-l, I ) +135 .*AUX IN -2 , T) + 

1 10 8.* AUX I N-3, I )+ AUX I N-4, I) > - . 0 23 437 5* I AUX I N+6, I ) + 1 8 . * AUX I N+5 , T > 
29.*AUX (N+4, I) )*H 
AUXIN-7, !)=AUX(N-o,T) 

631 AUXIN+4, I) = MJX(N + 5,n 
OE LT = X-H 
X= DEL T- ( H+H > 

T SW2 =9 
GOTH 700 

63? DO 637 1=1 , NOT M 
AUXIN-2, I 5 = Y( I ) 

AUX I N + 5, T)=DERY( I > 

637 Y I T ) = AUX I N-4 , T ) 
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X= X— ( H+H ) 

T SW2 = 1 0 
GOTH 700 

04 X =DE IT 

HO 63 6 I =1 ,NDI M 

OELT = MJX ( N+5, I ) + AUX(N + 4 t I) 

OELT = OELT4-OELT+nELT 

C AUV( 16, T) =8. 96 7960* ( AUXIN- l , I)-Y( T) ) — 0 . 7 6 1 1 11 *H* ( AUX( N+6, I ) +-OELT 
1+DERX (II) 

05 AUX(N+0,I)=DERY( I ) 

GOTO 606 

ENO 
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SUR ROUTINE GELG(R,A,M,N, EPS, TER) 
DIMFNSTON Ml ) ,R (1) 

TMM) 23,23, 1 

1 IFR=C 
PIV=0. 

M 

nm=n*m 
O n 9 I =i , mm 

TR=ARS( A ( L ) ) 

I F ( Tfl-P T V ) 3 , 3 , 2 

2 PTV=TR 
I = L 

3 CONTINUF 

tol=eps*ptv 

L ST= 1 

no 17 K= l , M 

IP(PI V) 23,23,4 

4 TF(IER)7,5,7 

5 IF (PIV-TPU 6, 4,7 

6 I F R= K — 1 

7 P T V I = 1 . / A ( I) 

J=( T-l ) /M 

1= T-J*M-K 
J= J+l -K 
DO R 1 =K,NM,M 
LL=L+T 

TR=P IV I*R ( LL ) 

R(LL)=R(L) 

R R ( L ) =TR 

IF(K-M)9, 18, 19 
9 LFND=LST+M-K 
I F ( J) 12, 12, IP 

10 I I = J *M 

DO 11 L=LST , LEND 
TR=A(L) 

LL = L + T I 
A(L)=A(IL) 

11 A ( LL ) =TB 

1? DO 13 L=LST , MM , M 
LL=L+T 

TB=P TVI*Ai LL ) 

A(LL )= A( L) 

13 A( L ) =TB 
A ( L ST ) =J 
PI V=0. 

LST=LST+1 

J=0 

DO 18 I I =LST, LFND 
P I V I =- A { I I ) 

I S T= I I +M 

J = J + 1 

DO 15 L=TST,MM,M 

LL =L- J 

A(L)=A(L)+PTVI*A(LL) 

TB=ARS( AIL ) ) 
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IF(TB-PIV) 15,15,14 
14 PTV=TR 
I =L 

1 5 CnNTTNUF 

DD 16 L =K» NM» R 
U = l + J 

16 P ( LL 1 =R ( LL ) ♦PIVI*R(L) 
1 7 1 ST=L ST+M 
1 R IF(M-1 >23,22,19 
19 TST=MM*.i^ 

L ST=M + 1 
DO 21 1=2, M 

I I=LST-T 
IST=I ST-LST 
L = T ST- M 
L=A(LH.5 
DO 21 J=IT,NM, m 
TR= R ( J > 

LL = J 

DO 20 K= I ST , MM, M 
LL=Ll ♦ 1 

2 n TR=TB-A( K )*R< LL ) 

K = J+l. 

R( J1=R(K) 

21 R ( K ) = TR 
2.7 RETURN 
?R I FR=- 1 
RETURN 
END 
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SUBROUTINE AFCT(X,A) 

THIS SUBROUTINE IS A PORTION OF PLPRV 

THESE DIMENSIONS MUST RF CHANGED IF THE INPUT 
POLYNOMIAL DIMENSIONS of °lpbv ape changed. 

DIMENSION P<3) ,0(3) 

MATRIX A MUST BE IN VECTOR FORM 

DIMENSION All) 

P,Q, AND N ARE PO|.Y NO M T AL COEFFICIENTS AND MUST BF 
MAIN PROGRAM AND PLACED IN COMMON AREA S . 

COMMON /S/P, Q» N 
C EPSILON IS ZERO VALUE FOR THIS SUBROUTINE 

FP S=1 .OF-07 
PI = P( l) 

P3 = Q ( 1 ) 

P 4=0 . 0 

IFIN-1 ) 10,11, 12 

11 P1=P1+P(2)*X 
P 4= Q ( ? ) 

GO TO 10 

12 P1=P1+P(2)*X 
P3=PB+?.0*P< 2 )+Q< 2 )*X 
P4 = Q I 2 ) 

X TMM=) .0 

no 20 1 = 2, N 

X IM=XTMM*X 
XT=XI M*X 
R I =FLO AT I I ) 

P1=P1+P ( 1+ 1 ) *XI 

P3=P3 + ?.0*RI*P(H-1)*XIM*-Q(I + 1)* : XI 

P4=P4 + RI*( (RI-l.O)APC t + 1 )*XTMM+Q< IM) *XIM) 

X I MM= X I MM* X 
20 CONTINUE 

in IF { ARS (PI J-EPSM B , IB, 14 
IB A ( 2) =0 .0 
A ( 3 ) = 0 . 0 
A ( 4 1 =0 .0 

TF(ABS(P3)-EPS)15,15,16 

15 A(1)=0.0 
RFTURN 

16 A(1)=-P4/PB 

return 

14 A(1)=0 .0 

A I 2) =-P4/Pl 
A ( B ) = l . 0 
A( 4) = -P3/Pl 
RETURN 
END 
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SUBROUTINE FCT ( X t F ) 

C THIS SUBROUTINE IS A PORTION OF PLPBV 

DIMENSION F(2) 

F ( 1 ) = n .0 
F { 2 ) = C . 0 
RETURN 
END 



SUBROUTINE OFCT(X,DF) 

C THIS SUPROUTINF IS A PORTION OF Pl.PBV 

DIMENSION DF ( 7 I 
DF ( 1 ) =0. 0 
DF( ?) =0.0 
RETURN 
END 



SUBROUTINE LOC ( I , J, T R, N, M, MS) 

IX=I 

JX=J 

IP(MS-l) 10 » 20 1 30 
10 TR X = N*( J X-l ) + T X 
C 7 O TO 36 

20 I F ( IX-JX) 2?t24»24 
2? trx=!X+( JX*JX-J X)/? 

GO TO 36 

2 A IP X= J X ■*■{ IX*IX-IX) /7 
GO TO 36 
30 IRX=0 

I F ( IX-JX) 3 6 » 32 1 36 
3? IRX = I X 
36 I R=I RX 
RETURN 
END 
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SUBROUTINE FCT { X , F ) 

C THIS SUBROUTINE IS A PORTION OF PLPRV 

DIMENSION F(2) 

F ( l ) = A . P 
F( 21=0.0 
RETURN 
END 



SUBROUTINE DFCT ( X , DF ) 

C TUTS SUBROUTINE IS A PORTION OF Pl.PBV 

DIMENSION DF ( 7 1 
DF ( l ) =0. 0 
DF ( 2 ) =0 . 0 
RETURN 
END 



SUBROUTINE LOC ( T , J, I R, N, M, MS) 

IX=I 

JX =J 

IP(MS-l) 10,20,30 
10 I R X= N* ( J X— 1 1 + T X 
GO TO 36 

20 I F ( I X— J X ) 22,24,24 
22 TRX=IX+( JX*JX-JX ) /? 

GO TO 36 

24 TP X= J X +( I X* I X— I X ) /2 
GO TO 36 
30 IRX=0 

IF(IX-JX) 36,32,36 
32 IRX = I X 
36 I R=T RX 
RETURN 
END 
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SUBROUTINE El GEM ( A , R , N , MV ) 

DIMENSION All ),R(1 » 

5 RANGE= 1. OE-6 

IF(MV-l) 10,25,10 
10 I 0=-N 

DO 20 J=1 » N 
1 0= IQ+N 
DO 20 1=1 , N 
I J=IQ + I 
R ( I J ) =0 . 0 
IF(I-J) 20,15,2° 

15 R ( IJ) = 1.0 
20 CONTI NUF 
25 ANOP M=0 .0 
DO 35 1 = 1, N 
DO 35 J=I,N 
IF(I-J) 30,35, 30 
30 IA=I+< J*J-J )/? 

ANORM= ANORM + AI I A ) *A( IA) 

35 CONTINUE 

IF(ANORM) 165, 155, AO 
AO AN0RM=I.A1A*SQRT( AN0RM1 

ANRMX = ANORM*RANGE/FLOAT( N) 

I ND = 0 

thr=anorm 

A5 THR=THR /FLOAT ( N) 

50 L = 1 
55 M=l+1 
60 MQ=(M*M-M)/2 
LQ=(L*L-L>/2 
(_ m=l + M Q 

62 I F ( ARS( A(LM ) 1-THR ) 130,65,65 

65 I ND=1 
LL=L +LQ 
mm=m+mq 

y=o.s*(A( LL)-A(MM) ) 

68 Y=-A(! M ) / SORT { A(l. M)*A(LM)+X*X) 

IF ( X ) 70,75,75 
70 Y=-Y 

75 S I NX = Y / SORT ( 2 • 0* ( 1.0+ ( SORT ( 1 ,0-Y*Y) ) ) ) 
SINX2=SINX*SINX 
78 C.OSX= SQRTI 1.0-SI NX 2) 

COSX2=COSX*COSX 
ST NO S =SINX*COSX 
I LQ=N* ( L-l ) 

IMQ=N* { M— 1 ) 

DO 125 1=1, N 
10= ( I* l-TI/2 
IF(I-L) 80,115,80 
80 if(I-M) 85,115,00 
fl<S IM=I + MQ 
GO TO 95 
90 I M=M+ I Q 

95 IF(I-L) 100,105,105 
ion IL=I+LQ 
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GO TO ll'' 

105 I L=L+I 0 

110 X=A(IL)*COSX-A( I M) *S INX 

A(IM)=A(IL)*SINX+A( T M ) *COS X 
AIR) =X 

115 TF(MV-i) 120,125,120 
12 n ILR= ILQ+ I 
I MR= I MQ+T 

X=R( ILR)*COSX-R( IMR)*SINX 
R (IMR) =R( TLR)*5 INX+Rl IMP ) *COSX 
R ( ILR )=X 
125 CONTINUF 

X=2.0*A(LM) *SINCS 

Y=A( Ll )*C0SX2+A( MM) *$I NX2-X 

X=A(LL )*SINX2+A(MM)*C0SX?+X 

A(LM)=(A( LL)-A(MM) )*STNCS+A(LM)*(C0SX2-SINX2 ) 
A(LL)=v 
A ( MM ) -=X 

1^ IF(M-N) 135,140,155 
13 5 M= M+ 1 

GO TO 60 

140 IF(L-(N-1)) 145,150,145 

145 L=L+ 1 

GO TO 55 

15° I F ( TNO-1 ) 160, 155,160 
155 TND=0 

GO TO 50 

160 IF(THR-ANRMX) 165,165,45 
165 IQ=-N 

On 1R5 T=1,N 

I Q=I Q+ N 

LL=l+( 1*1-1 ) /? 

JQ=N*C 1-2) 

DO 18 5 J=I,N 
J 0= J Q+ N 

M M= J + ( J* J — J ) /2 
IFIAILL)-A(MM) ) 170,135,155 

170 X = A 1 LL ) 

A(LL )=A(MM) 

A ( MM ) = X 

IF(MV-l) 175,135,175 
175 DO 1 80 K=1 , N 
ILR= TO+K 
I mr= jo+K 
X =R ( ILR ) 

R ( ILR ) =R ( IMP ) 

180 R ( IMR )=X 
185 CONTINUE 
RFTUPN 
END 
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APPENDIX A6 



SUBROUTINE PLPBV 

This appendix presents a complete listing of a suggested format 

and some of the computation statements necessary to solve the PLPBV 

model which is a special case of the LPCM, As stated in the subroutine 
description, the program presented here is incomplete and represents 
at best a format for further use developing the complete solution to 
the LPCM utilizing the integral equation methods of this thesis. 

One important result of the use of this subroutine is that the 
total computation time necessary for solution is very short. Sub- 
routine PLPBV, the MAIN program calling it, and all of the subroutines 
used compiled under WATFOR in 12 seconds and executed completely in 
15 seconds. More elaborate programs would undoubtedly require more 
than 15 seconds computation time, but the fact that a problem as general 

as PLPBV can be solved in a time this short is very encouraging for 

more general versions of the LPCM. A flowchart of PLPBV is presented 
in Chapter L2. 

For completeness, the MAIN program used to call PLPBV as a test 
is presented at the end of this appendix. In general, this main 
program would employ the approximations and combinations of upper and 
lower solutions for a complete distillation column modeling using 
PLPBV to solve the upper and lower versions of the LPCM, 
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SUBROUTINE PL PBV I PRMT , P , Q , N, BC,F,US ,UT,U ) 



C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



NOTE ***#*«**************»>**********#***************** 
THIS PROGRAM IS INCOMPLETE IN SEVERAL ASPECTS AND 
REQUIRES EXTENSIVE WORK TO BE USEFUL IN GENERAL 
PURPOSE COMPUTATION. 

SUBROUTINE PLPBV IS PRESENTED HERE AS A FORMAT FOR 
AND EXAMPLE OF THE STEPS NECESSARY TO APPLY THE 
INTEGRAL EQUATION TECHNIQUE OF THIS THESIS. 

NOTE #*$*##**<“!'*****$*$*******£***#****+*****#******** 



**=#*#* ijtJj!####*#**#.##*#*###*#**#**:*#*#**#**##### ** $*«$!(£* 



SUBROUTINE PLPBV < PRMT , P , Q, N , BC , F , US , UT ,U > 

PURPOSE 

TO SOLVE A LINEAR, SECOND ORDER, P AR ABOL I C, 
PARTIAL DIFFERENTIAL EQUATION WITH BOUNDARY 
CONDITIONS AND POLYNOMIAL COEFFICIENTS, 
D/DXX{P(X)U(X»T))+D/DX(Qfx)UIX»T) ) =D/DT I U( X , T ) ) 
WHERE P AND Q ARE POLYNOMIALS IN X, IN RESPONSE 
TO STEP INPUTS, WHERE F IS STEADY STATE AT T=TL. 



USAGE 

CALL 



PL PBV IPRMT ,P ,Q,N,BC,F,US,UT,U) 



DESCRIPTION 

F 

BC 



PRMT 



PRMT ( I > — 
PRMT 11)- 
PRMT (3)- 
PRMT 14)- 
PRMTI5 )- 
PRMT {&)- 
PRMT { 7 ) - 
PRMT (8)- 
PRMTC9)- 
PRMT 1 105 



OF PARAMETERS 

INPUT STEADY-STATE DISTRIBUTION AT TL 
INPUT BOUNDARY CONDITIONS OF THE FORM 
3C ( 7 ) U+ BC ( 9 ) UX=BC ( 1) AT X=XU 
3C<4)U+BC<6) UX=BC( 2) AT X=XL 
BC MUST BE AT LEAST OF DIMENSION 10 
PLPBV SETS BCI 3, 5,8,101=0 FOR LBVP USE 
INPUT VECTOR WHICH SPECIFIES THE PARA- 
METERS OF THE X AND T INTERVALS AND OF 
THE ACCURACY FOR SUBROUTINES USED BY 
PLPBV, MUST BE AT LEAST DIMENSION 10 

THE X VARIABLE 
THE X VARIABLE 
OF THE X VARIABLE 
FOR X 



LOWER BOUND XL OF 
UPPER BOUND XU OF 
SPATIAL INCREMENT 
UPPER ERROR BOUND 
TERMINATION PARAMETER 
LOWER BOUND TL OF THE 
UPPER BOUND TU OF THE 
TIME INCREMENT OF THE 
ERROR WEIGHT LTE 1.0 



FOR X 

T VARIABLE 
T VARIABLE 
T VARIABLE 
FOR LBVP USE 



ERROR WEIGHT LTE 1.0 FOR LBVP USE 
BOTH ERROR WEIGHTS ARE USUALLY 1.0 
INPUT VECTOR SPECIFYING THE 
COEFFICIENTS OF THE FIRST 
POLYNOMIAL OF DEGREE N 
INPUT VECTOR SPECIFYING THE 
COEFFICIENTS OF THE SECOND 
POLYNOMIAL OF DEGREE N 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



c 

c 

c 

c 

c 

c 

c 



N - INPUT SPECIFYING THE MAXIMUM 

DEGREE OF P AND Q 
THE COEFFICIENTS MUST BE IN THE 
ORDER OF LOW TO HIGH INDEX 
CORRESPONDING TO LOW TO HIGH 
DEGREE IN EACH POLYNOMIAL 
BOTH P AND 0 MUST BE AT LEAST OF 
DIMENSION N WITH THE REMAINING 
COEFFICIENTS IN EITHER P OR Q SET=0 

US - OUTPUT STEADY-STATE SPATIAL 

DISTRIBUTION VECTOR 

UT - OUTPUT TRANSIENT RESPONSE 

U - OUTPUT U=US+UT OVERALL SOLUTION TO THE 

TWO POINT BOUNDARY VALUE, STEP INITIAL 
VALUE PROBLEM. 

REMARKS 

(1) STEADY-STATE PORTION WRITTEN FOR MAX N=3, 

THIS CAN BE EASILY CHANGED BY USING 
DIFFERENT DIMENSION STATEMENTS ON P AND Q 

(2) TRANSIENT PORTION WRITTEN FOR P(1),Q(1), 

AND 0(2), CHANGING THIS REQUIRES ANALYTICAL 
APPLICATION OF THE INTEGRAL EQUATION 
TECHNIQUE PRESENTED IN THIS THESIS. 

(3) THIS SUBROUTINE HAS BEEN INITIALLY WRITTEN 
TO EVALUATE AT TEN POINTS IN X AND AT NINE 
POINTS IN T. THIS COULD BE CHANGED TO A 
MORE GENERAL METHOD IF DESIRED. 

(4) THERE ARE SEVERAL ERROR OUTPUTS IN PLPBV. 

15) SUBROUTINE PLP8V HAS NOT BEEN OPTIMIZED. 

(6) SUBROUTINE PLPBV HAS BEEN WRITTEN MAINLY TO 

DEMONSTRATE THE STEPS NECESSARY TO APPLY 
THE INTEGRAL EQUATION SOLUTION OF THIS THESIS. 

17) IT IS EXPECTED THAT SEVERAL 'BUGS* REMAIN 
IN SUBROUTINE PLPBV AS WRITTEN HERE. 

( 8 } AS THE AFCT OF APPENDIX A5 IS NOW WRITTEN, 

A STATEMENT, COMMON /S/P,Q,N ,IS REQUIRED 
IN THE MAIN PROGRAM. THIS COULD BE CHANGED 
BY INCREASING THE PARAMETER DIMENSIONS IN 
SUBROUTINE AFCT. 

SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 

SUBROUTINE PLPBV REQUIRES IBM/SSP SUBROUTINES- 
LBVP, GELG » EIGEN, AND LOC 

AND USER FURNISHED SUBROUTINES- 

AFCT , DECT, FCT , AND OUTP FOR USE IN LBVP 

# *****5}C# **##***# «#*##**#=('*#*♦*# **♦«*#*** 

DIMENSION PRMT (1) ,P( 1) ,Q(1) , BCIIItFIl ) ,US( 1) ,UT( 1) ,U( II 
DIMENSION D ( 4 ) , AK ( 100 ) , UZ I 1C I , X II 1 00) , XJ( 1 00 ) 

DIMENSION AKSI55) , THC 10 > , E IG 1 100 ) , TL AM 1 10 ) ,W II 100) 



INPUT AND PARAMETER SETUP SECTION 
ND I M=2 
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NV= 10 
NT=9 
CN= 10 .0 
C T=8 . 0 
E PS= 1 . 0E— 05 
A A= PRMT ( 1 ) 

BB=PRMT ( 2 ) 

TL=PRMT (6) 

THH=PRMT (7) 

I FI ABS t P ( 1) )-EPS >33,33,34 

33 WR ITE ( 6*35 ) P( 1 ) 

35 FORMAT I IX,' p{ 1 ) INVALID, P(l> = *,E15.7,//> 

RETURN 

34 SPZ = SQRT { P( 1 ) ) 

SSPZ=SORT(SPZ ) 

STEADY STATE SECTION, USES SUBROUTINE LBVP 
STEADY-STATE VALUES RETURNED THROUGH XPAR 
DIMENSION Y ( 2 J , DERY ( 2 ) , AV { 4 ) , AUX C 20 , 2 ) 

DIMENSION RI2) ,8(4} ,C(4),TPAR( 15) , XPAR (15) 

EXTERNAL FCT , DECT, AFC T, OUTP 
DO 55 1=1,5 
XPAR { I )=PRMT ( I ) 

TPARI I )= PRMT { 1 + 5 ) 

55 CONTINUE 
R(1)=BC(1) 

R { 2 ) =BC( 2 ) 

BC (3 ) =0.0 
BC ( 5 ) =0 .0 
BC ( 8 )=0.0 
BC { 10) =0.0 
DC 56 1 = 1,4 
B( I ) =BC ( 1+2 ) 

C ( I ) = BC1 1+6) 

56 CONTINUE 
DERYI1)=PRMT(9> 

DERY i 2 )=PRMT 110) 

CALL LBVPJXPAR , 8 , C, R, Y , DE RY, NDI M , I HL F , AFCT , FCT , DECT , 

1 OUTP , AUX, AV) 

TESTS APPLIED UPON RETURN FROM LBVP 
IF(IHLF-13) 39,40,41 

41 WRITEI6.42) 

42 FORMAT! IX, » LBVP HAS IHLF=14, NO SOLUTION ',//) 

RETURN 

39 I F C I HLF- 11)44,45, 40 

45 WR I TE ! 6 » 46 ) 

46 FORMAT { IX, ' LBVP HAS I HLF GT 10, NO SOLUTION',//) 

RETURN 

40 WR IT E ( 6 * 43 ) 

43 FORMAT ( IX ,' LBVP HAS IHLF=13 OR 12, PARAMETER ERROR',//) 
RETURN 

44 WRITE(6,47) I HLF 

47 FORMAT (IX,* LBVP RETURN SATISFACTORY, IHLF= *,I3,//) 



C 
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TRANSIENT SECTION, USES SUBROUTINES EIGEN AND GEL G 
TRANSFORMING BOUNDARY CONDITIONS 

D( 1) = BC(7)— ( BC(9)*(Q( 1 ) +0 { 2 J *PRMT ( 2 ) ) )/(2.0*P( 1)) 
D(2)=BC(4)-(BC(6)*(Q(1>+Q(2)*PRMT( 1) ) )/(2.0*P(l>) 

D( 3 ) — BC ( 9 )/S PZ 
D ( 4 ) = BC ( 6 ) /SPZ 
DTEST-=D( 1 )*D(4)-D(2)*D(3) 

C0NE=D(1)*D(2)— OTEST 
IFCABSICONE )— EPS )31»31,30 

31 WRITE (6,32) CONE 

32 FCPMAT (IX, • Cl IS • , E 1 5 . 7 , • SE TT I NG Cl TO l.E-05 ',//) 
CONE= EPS 

SET UP INTERVAL FOR KERNEL EVALUATION 
30 DELTX=(PRMT(2)-PRMT(1) )/CN 
CC=DELTX*CN/P( 1) 

PROGRAM RUNS FOR 10 OELTX AND 10 DELTZ 
FOR GREATER ACCURACY MORE INCREMENTS COULD BE USED 
DELTZ=CC /CN 

BEGIN KERNEL EVALUATIONS 
EMZ=1. 5*0(2) +0.2 5*Q( 1)*Q(1)/P( 1) 

C 

C SEVERAL ALTERNATIVE KERNEL EVALUATION ROUTINES 

C MUST BE DESIGNED AND PLACED IN THIS SECTION FOR 

C CASES SUCH AS ** 

C (1) M ( 0 ) =0 

C (2) M(Z) = M(0) 

C ANC OTHERS. 

C THESE ROUTINES MUST USE EQUATIONS A1.23 AND A1.24 

C OF THIS THESIS. 

C 

C TEST M (0 ) 

IF(A8S(EMZ)-EPS)50»50, 51 

50 WRITE(6,52> EMZ 

52 FORMAT ( 1 X, * M(0) =0 , INVALID FOR PLPBV, M ( 0 ) = ' , E 1 5 . 7 , / / ) 
RETURN 

51 DK=D( 1 ) +D ( 3 ) 

FP=4.0*P( 1 ) 

LK=0 

L S=0 

DO 49 JK=1 , NV 
DO 48 I K=1 , NV 
C COMPUTE K4(Z,S) 

LK=LK+1 
Z= I K*OEL T Z 
S=JK*DELTZ 
XS= SPZ*S + AA 
CMS=CC-S 

” DMD’Z= C D'( 4 J — D ( 2 ) *Z ) /CONE 
XCMS=SPZ*CMS+ AA 

EMS=EMZ+(2.0*0(1)*0(2)*XS+0(2)*0(2)*XS*XS)/FP 
XQ=G( 2 ) *XCM S 
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EMCMS=EMZ+(2.0*Q ( 1 ) *XQ+XQ*XQ ) /FP 
KF=DMDZ*(D( 1)*EMCMS+D(3)*EMS-DK*EMZ)/EMZ 
IF ( IK-JK) 54,54,53 
54 AK ( L K } =KF 
LS = LS+1 

AKS!LS)=KF Note: obvious error 

GO TO 48 KF should be RKF, real. 

53 ZMS=Z-S 

XZMS=SPZ*ZMS+AA 
X0Z=0 12 ) *XZMS 

AK ( LK ) = ( (2, 0*Q! 1 )*XQZ+XQZ*XQZ )/< EMZ*FP) )+KF 

48 CONTINUE 

49 CONTINUE 
C 

C COMPUTE EIGENVALUES AND EIGENVECTORS OF AKS 

M V=0 

CALL EIGEN! AKS, E IG, NV , MV ) 

TAKE EIGENVALUES TH(I) FROM DIAGONAL OF AKS { I ) 

EPTH= 1 . OE— 10 
I S= 1 

DO 60 1=1, NV 
TH( I )=AKS< IS) 

I S=IS + ( 1 + 1) 

IF(TH( I )— EPTH)57,57,58 

58 TLAM(I)=1.0/(TH( I)#DELTZ> 

GO TO 60 

57 WRITE(6,74)TH{ I), I 

74 FORMAT! IX, ' THETA! I ) INVALID, T H( I ) = • , El 5 . 7, • 1= • , I 3> 

WRITE!6,59) 

59 FORMAT! IX, • SETTING TH! I )=EPTH* ,/) 

TL AM! I)=EPTH 

60 CONTINUE 

THIS PROGRAM HAS 8EEN SET UP SO THAT DELTX AND DELTZ 
CORRESPOND. THIS MAKES W I ! Z I ) =W I ( X I ) . IF THIS WERE 
NOT DESIRED, THEN A SECTION UTILIZING EQUATION A3. 10 
WOULD HAVE TO BE EMPLOYED. 

COMPUTE WIIX) FROM WI(Z), THETA! I )=TH (1) , WI ! Z ) = E IG! Z ) 
IJ = 0 

DO 63 1=1, NV 
DO 64 J= 1 , N V 
IJ=IJ+1 

W I! I J )=EIGi I J) 

64 CONTINUE 
63 CONTINUE 

TRANSFORM WI (X) TO XI! X) 

X=DELTX+AA 
DO 61 J=1 ,N V 
XMA=X— AA 
XPA=X+AA 

XFUN=( XMA*Q(1) )/(2.0*P( 1) ) + ! XMA*XPA*Q ( 2 ) ) /FP 
EXX=EXP!-XFUN> *SSPZ 
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X=X+DELTX 
DO 62 I = 1 » N V 
JI=J+NV*( 1-1 J 
XI( JI )=EXX*WI( JI ) 

62 CONTINUE 
61 CONTINUE 
C 

C BEGIN EIGENFUNCTION EXPANSION TO MEET INITIAL CONDITION. 

DELTT=(THH-TL) /CT 
T = TL 
C 

C FORM STEADY-STATE VECTOR 

DO 65 1=1, NV 
UZ{ I ) =F ( I )-US ( I ) 

65 CONTINUE 
C 

C FORM EXPANSION M ATR I X , X J ( X ) , AT T=TL 

IJ=0 

DO 66 1=1, NV 

ext=exp(~tlam( i)*t> 

DO 67 J= 1 ,N V 
IJ=IJ+1 

XJ( IJ)=XI( I J )#£XT 

67 CONTINUE 

66 CONTINUE 
C 

C SOLVE SIMULTANEOUS EQUATIONS FOR EIGENFUNCTION 

C CONSTANTS USING SUBROUTINE GELG 

NG= 1 

CALL GELG(UZ,XJ.NV,NG,EPS, I E R > 

C 

C TEST IER UPON RETURN FROM GELG 

IF< IER )36»37,38 

36 WRITE<6,68> 

68 FORM AT { 1 X , * NO RESULT, PIVOT ELEMENT=0 IN GELG *,/) 

RETURN 

38 WR I TE 1 6 » 69 ) 

69 FORMATI1X, • POSSIBLE LOSS OF SIGNIFICANCE IN GELG*,/) 

C 

C THIS SECTION SETS UP U(X,T) IN A FORM ACCEPTIBLE 

C FOR PLOTTING BY SUBROUTINE PLOT OF APPENDIX A4. 

C THE FIRST COLUMN Or U(I) IS THE INDEPENDENT 

C VARIABLE X, THE SECOND COLUMN IS THE STEADY-STATE, 

C AND THE REMAINING 8 COLUMNS ARE RESPONSES. 

37 X= A A 

DC 70 1=1, NV 
X=X+DELTX 
UU ) = X 

U< 1 + 1 OJ = F { I ) 

70 CONTINUE 
C 

C CALCULATE TRANSIENT SOLUTION MATRIX UTCX.T) 

C 

C SET TIME-SPACE GRID 

I JK=20 
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T=7L 

NTC=NT— 1 
DC 71 1=1, NTC 
T= T +DELTT 
X=AA 

DO 72 J=1,NV 
X=X+DELTX 

TIME-SPACE GRID NOW SET 



73 

72 

71 



SUM SERIES FOR EACH VALUE 
I JK= I JK+1 
UT < I JK ) = 0.0 
DO 7 3 K=1,NV 
I X= J+NV* ( K— 1 ) 

EXT=EXP(-TLAM{ K)*T) 

UT( I JK)=UZJK)*XI ( IX ) *EXT +UT ( I JK ) 
CONTINUE 

U( IJK >=US( JJ+UTi I JK) 

CONTINUE 

CONTINUE 

RETURN 

END 
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HAYES*KP=29*TIME=2«PAGES=50 



S JOB 

MAIN PROGRAM TO USE PLPBV AS A TEST 
DIMENSION BC{ 10) ,Ut 100),US( 10) ,UT( 100) ,F{ 10) 

DIMENSION P C 3 ) *0(3) *A( 2*4), PRMT (20) 

COMMON /S/P , Q , M 

PLPBV INPUT SECTION 

M= 3 

READ (5, 100 ) (PI I ) * Q! I) , 1 = 1 * M ) 

100 FORM AT ( 2 F10 . 5 ) 

READ! 5» 101) l(A(I,J),J=l,4),I=l,2) 

101 FORMAT ( 4F1G .5 ) 

READ (5.102) (PRMT( I) ,1=1, 10) 

102 FORMAT (2E15. 7) 

BC ( 7) =A{ 1,1) 

BC (9 ) =A{ 1 ,2 ) 

BCI 1)=A( 1,4) 

BC(4)=A(2,1) 

BC(6)=A(2,2) 

BC ( 2 ) = A ( 2,4) 

SET UP LINEAR STEADY-STATE AS A TEST 
X = 0.0 

DO 10 KSS = 1 , 10 
F l KSS ) =G . 0 
X=X+0 . 1 
US(KSS) = X 
10 CONTINUE 
C 

C TEST STEADY-STATE INPUT FOR PLPBV FROM UC1, APPENDIX A4 

PRMT ( 1 1 ) =0. 06 
PRMT ( 12)=0. 126 
PRMT ( 13) =0. 200 
PRMTC 14 ) =0 • 27 
PRMT( 15 ) =0 . 37 
PRMT! 1 6) =0. 47 
PRMT { 17 ) =0 » 59 
PRMT! 1 8 ) =0, 7 1 
PRMT! 19 ) =0. 84 
PRMT I 20 ) = 1 • 00 

CALL PLPBV! PRMT ,P ,Q, M, BC , F, US, UT,U ) 

CALL EXIT 
END 
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SECTION 6 



BIBLIOGRAPHY AND REFERENCING (B) 

B1 BIBLIOGRAPHY 

B2 BIBLIOGRAPHY REFERENCES LISTED BY AREA OF APPLICATION 
B3 BIOGRAPHICAL NOTE 



AN UNDERSTATEMENT: 

" THE LITERATURE ON DISTILLATION IS 

R. 



VOLUMINOUS - " 

J. HENGSTEBSCK (H-18) 



ONE OF THE MAJOR EFFORTS OF THIS THESIS TURNED OUT TO BE THE 
COMPILATION OF THE BIBLIOGRAPHY PRESENTED IN THIS SECTION. THE 
REFERENCES LISTED APPLY TO FOUR MAJOR AREAS PERTINENT TO THIS 
THESIS. 

1. GENERAL THEORY OF DISTILLATION 

2. DISTILLATION COLUMN DYNAMICS 

3. DISTILLATION COLUMN CONTROL 

4. MATHEMATICS AND COMPUTATION 

THE AUTHOR HOPES THAT UTILIZATION OF CHAPTER B2, BIBLIOGRAPHY 
REFERENCES LISTED BY AREA OF APPLICATION, WILL RESULT IN CONSIDER- 
ABLE SAVINGS OF SEARCHING TIME AND EFFORT FOR ANYONE INTERESTED IN 
THESE AREAS. 



- 234 - 



CHAPTER B1 



BIBLIOGRAPHY 

This chapter presents a bibliography of 352 references pertinent 
to the general area of distillation and to the specific areas of this 
thesis. The references consist only of those in English, although 
extensive literature on distillation has been published in the foreign 
journals, especially in Russian and German. For the most part, the 
references were taJcen from the following journals for the years from 
about 1955 to 1969. 

1. Industrial and Engineering Chemistry 

2. American Institute of Chemical Engineers 

3. Chemical Engineering Progress (and Symposium Series) 

4. Transactions of the Institution of Chemical Engineers 

5. Chemical Engineering Science 

6. British Chemical Engineering 

7. Canadian Journal of Chemical Engineering 

The entries in this bibliography are in the order of the first 
letter of the author's last name, for the reader's mnemonic convenience, 
but no attempt has been made to subalphabetize within each group for 
this author's convenience. 

It is this author's intention that Chapter B2 - Bibliography 
References Listed by Area of Application be used in conjunction with 
this chapter for any given subject or area of research. In each ref- 
erence in this bibliography an attempt has been made, especially with 
journal articles, to present as complete a description as possible of 
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the reference, even at the cost of complete overspecification and 
possible redundancy of information. 

The criterion of availability of each reference in the M.I.T. 
libraries (except for theses) placed a significant constraint on the 
number of references which have been listed here. The author chose 
this as his stopping point in the compilation of this bibliography. 
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CHAPTER B2 



BIBLIOGRAPHY REFERENCES LISTED BY AREA OF APPLICATION 

Ihe purpose of this chapter is to classify the references in 
Chapter Bl, the Bibliography, according to their principal areas of 
application. The general format of this classification is presented 
in Table B2.1. Many of the references apply to more than one area and 
are so listed for the user's convenience. Rather than merely list the 
references under each area in alphabetical order, the author has chosen 
to present them column by column in the order of decreasing utility to 
the study of the subject or of decreasing clarity. In other words, the 
first listed in each area seem to the author to be the most important 
for anyone researching that area. Recognizing that such a listing is, 
indeed, very subjective, the author apologizes to anyone who may find 
them "out of order" with respect to his particular slant on the subject. 
B2.1 General Theory of Distillation 
Textbooks 

Extensive Bibliographies and Literature Surveys 
References of Historical Interest 
General Distillation 
Dynamic or Transient Analyses 
Steady-State Analysis and McCabe-Thiele Diagrams 
Structural Design 
Economics and Operations Analysis 
Thermodynamics 
Hydrodynamics 
Chemistry 

Philosophy Table B2.1 (Contd. ) 



B2.2 Distillation Column Dynamics 



Textbooks 

Iheses 

Reviews, Bibliographies, and Literature Surveys 
Dynamic Models or Solutions 
Discrete Plate Equations 

Frequency Analysis or Laplace Transform Solution 
Transient or Time Analysis 
Numerical Solution 
Digital Computer 
Hybrid Computer 
Analog Computer 
Analytical Analysis 
Continuous Spatial Equations 

Frequency Analyses or Laplace Transform Solution 
Transient or Time Analysis 
Experimental Transient Behavior 
Frequency Response 
Time Response 
Cyclic Distillation 
B2, 3 Distillation Column Control 
Textbooks 
Iheses 

Extensive Bibliographies and Literature Surveys 
Conventional Control Systems 
Digital Control 

Hybrid Control __ , . , , . 

Tabl e B2 . 1 ( Contd. ) 
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Analog Control 
Instrumentation 



Control Systems Using Dynamic Models 
Digital Computer Control 
Hybrid Computer Control 
Analog Computer Control 
Optimal Control 
Distributed or Modal Control 
B2.4 Mathematics and Computation 

Ordinary Differential Equation Theory 

Partial Differential Equation Theory 

Integral. Equation Theory 

Mathematical Transformations 

Matrix Mathematics 

Numerical Solution Techniques 

Boundary Value Problems 

Eigen - Values, Vectors, and Functions 

Special Functions 

Computation and Computer Programming 
Table B2.1 - FORMAT FOR CLASSIFICATION BY AREA OF APPLICATION 
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B2.1 GENERAL THEORY OF DISTILLATION 



Textbooks 



(B-l) 


(A-19) 


(H-18) 


(R-22) 


(T-l) 


(G-3) 


(B-15) 


(L-26) 


(R-17) 


( B-12 ) 


(H-10) 


(M-14) 


(C-5) 


(S-8) 


(H-5) 


(R-23) 


(0-2) 


(C-l) 


(S-9) 


(C-7) 


(R-21) 


(H-2) 


(P-3) 


(s-3) 


(M-2) 


(V-2) 


(H-9) 


(P-10) 


(S-l) 


(V-3) 


Extensive 


Bibliographies and Literature Surveys 


(W-l6) 


(G-9) 


(w-13) 


(W-18) 


(P-6) 


(B-17) 


(R-12) 


(w-15) 


(W-26) 


(G-14) 


(B-18) 


(P-11) 


(P-7) 


(2-3) 


(W-23) 


(B-19) 


(R-8) 


(R-25) 


(W-8) 




(W-10) 


(H-7) 


(R-18) 


(W-17) 




(w-ii) 


(W-9) 


(H-8) 


(G-10) 




(P-8) 


(W-12) 


(R-17) 


(w-19) 





References of Historical Interest 
General Distillation 



(L-19) 


(L-25) 


(V-2) 


(P-3) 


(P-12) 


(G-14) 


(A-17) 


(M-10) 


(M-9) 


(R-29) 


(C-15) 


(T-6) 


(1-1) 


(u-i) 


Dynamic or 


Transient Analysis 








(M-8) 


(B-29) 


(L-2) 


(B-7) 


(R-2) 


(R-12) 


(J-1) (F-5) 



Steady State Analysis and McCabe - Thiele Diagrams 



(B-l) 


(H-22) 


(C-l) 


(R-21) 


(S-23) 


(L-25) 


(B-7) 


(R-18) 


(M-10) 


(B-10) 


(P-3) 


(H-30) 


(S-12) 


(R-30) 


(S-29) 


(A-19) 


(E-4) 


(M-18) 


(F-4) 


(J-7) 


(C-14) 


(T-l) 


(W-30) 


(P-12) 


(H-10) 


(P-8) 


(V-2) 


(L-12) 


(D-2) 


(E-2) 


(P-9) 


(T-7) 
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B2.2 



Steady State Analysis and 


McCabe - 


Thiele 


Diagrams (Contd.) 




(H-18) 


(V-2) 


(F-15) 


(A-10) 


(0-2) 


(s-15) 


(R-l) 


(R-l) 


(H-26) 


(s-15) 


(G-3) 


(A-17) 


(T-6) 


(s-5) 


(H-21) 


(R— 29) 


(C-9) 


(S-20) 


(H-5) 


(B-24) 


(P-13) 


(s-3) 


(I-D 


(s-l6) 


(M-8) 


(Z-2) 


(H-2) 


(B-26) 


(R-5) 


(S-24) 


(L-24) 


(s-1) 


(M-9) 


(R-26) 














Structural 


Design 














(L-13) 


(R-22) 


(F-10) 


(J-3) 


(M-ll) 


(L-12) 


( D-ll ) 




Economics and Operations Analysis 










(L-15) 


(M-13) 


(B-15) 


(S-8) 


(F-9) 


(F-10) 


(C-7) 


(v-3) 


Thermodynamics 














(H-9) 


(R-23) 














Hydrodynamics or Fluid Mechanics 










(F-4) 


(V-2) 


(R-23) 


(S-5) 


(P-9) 


(B-4) 


(S-19) 




(L-26) 


(B-12) 


(H-14) 


(G-U) 


(F-5) 


(M-24) 


(T-7) 




Chemistry 
















(B-16) 


(B-l) 


(H-17) 


(T-6) 


(s-5) 


(R-23) 






(H-8) 


(H-10) 


(F-9) 


(C-8) 


(H-20) 


(T-6) 






Philosophy 
















(A-13) 


(B-28) 














DISTILLATION COLUMN DYNAMICS 










Textbooks 
















(M-l) 


(H-7) 


(G-3) 


(F-16) 










Theses 
















(F-l4) 


(W-21) 


(D-14) 


(C-10) 


(B-32) 


(W-22) 


(S-27) 




(M-15) 


(M-12) 


(M-17) 


(Q-l) 


(0-3) 


(G-5) 


(W-20) 




(S-28) 


(A-5) 


(r— 6) 


(R-27) 


(S-22) 


(M-19) 


(W-24) 




(S-17) 


(C-4) 


(A -26^ 


(L-17) 


(S-18) 


(M-24) 







Reviews, Bibliographies, and Literature Surveys 



(A-3) 


(H-7) 


(R-25) 


(W-15) 


(z-3) 


(s-33) 


(L-3) 


(R-12) 


(B-17) 


(R-8) 


(W-12) 


(G-9) 


(w-13) 


(T-2) 


(W-l6) 


(B-18) 


(W-10) 


(W-ll) 


(B-19) 


(W-24) 


(W-23) 



Dynamic Models or Solutions 
Discrete Plate Equations 



Frequency Analysis or Laplace Transform Solution 



(A-3) 


(H-7) 


(R<25) 


(w- 15 ) 


(Z-3) 


(S-33) 


(L-3) 


(R-12) 


(B-17) 


(R-8) 


(w- 12 ) 


(G-9) 


(W-13) 


(T-2) 


(W-l6) 


(B-18) 


(w- 10 ) 


(w- 11 ) 


(B-19) 


(W-24 


(W-23) 


Transient 


or Time 


Analysis 









Numerical Solution 
Digital Computer 



(D-4) 


(R-16) 


(P-2) 


(T-9) 


(S-16) 


(B-32) 


(R-9) 


(P-17) 


(w- 25 ) 


(R-19) 


(R-2) 


(A-4) 


(D-9) 


(Y-l) 


(L-27) 


(D-8) 


(M-3) 


(L-18) 


(P-5) 


(S-20) 


(H-3) 


(R- 10 ) 


(D-6) 


(G-4) 


(W-4) 


(L-5) 


(R-ll) 


(A-9) 


(R-5) 


(R-l8) 


(R-20) 


(T-2) 


(s-6) 


(D-ll) 




(D-15) 


(w-i) 


(s-33) 


(R-l) 




Hybrid Computer 








(F-l) 


(F-l 3) 


(R-6) 


(P-12) 




Analog Computer 








(B-23) 


(P-D 


(H-3) 


(G-12) 


(L-5) 
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Analytical Analysis 



(B-3) 


(M-l) 


(S-31) 


(B-27) 


(H-20) 


(G-2) 


(R-7) 


(W-2) 


(G-2) 


(H-21) 


(A-3) 


(R— 8) 


(F-16) 


(z-3) 


(B— 7 ) 


(M-8) 


(R-32) 


(G-3) 


(P-5) 


(B-12) 



Continuous Spatial Equations 



Frequency Analysis or Laplace Transform Solution 



(J-6) 


(M-27) 


(H-4) 


(S-7) 


(w-5) 


(D-l) 




Transient 


or Time 


Analysis 










(j-i) 


(R-8) 


(M-l) 


(L-26) 


(T-5) 


(J-5) 


(W-14) 


(H-4) 


(D-3) 


(M-28) 


(P-7) 


(G-3) 


(L-20) 


(B-25) 


(K-2) 


(S-13) 


(K-3) 


(R-14) 


(H-25) 


(B-20) 


(C— 8 ) 


(0-1) 


(D-14) 


(G-2) 


(B-29) 


(J-2) 


(B-24) 





Experimental Transient Behavior 
Frequency Response 

(H-15) (H-23) (A-6) (H-24) (W-28) 

Time Response 

(H-3) (L-3) (B-3) (B-3) (H-15) (S-4) (D-12) (H-30) 

Cyclic Distillation 

(M-19) (A-14) (M-20) (B-17) (B-18) (S-4) 

B2.3 DISTILLATION COLUMN CONTROL 



Textbooks 

(A-12) (G-3) (B-13) (K-5) (C-l) (C-15) (A-ll) (L-l) 

Theses 

(B-2) (B-32) (M-5) (G-l) (S-18) (M-25) 
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Extensive Bibliographies and Literature Surveys 



(A-3) (B-17) (R-25) (R-12) (W-10) (W-20) (W-12) (W-2l) 

(W- 16 ) (B-18) (R- 8 ) (S-33) (w-ll) (b-19) (W-15) (V-23) 

Conventional Control Systems 
Digital Control 

(S-32) (H-29) (B- 6 ) (L-15) (L-24) (T- 8 ) 

Hybrid Control 



(P-13) (P-1) (F-12) 

Analog Control and Instrumentation 



(H- 16 ) 


(L-16) 


(B-10) 


(C- 15 ) 


(G-12) 


(s-30) 


(M-29) 


(S-9) 


(R- 31 ) 


(B-13) 


(B-31) 


(P-6) 


(K-4) 


(C-12) 


(P-14) 


(w-23) 


(G-3) 


(B-22) 


(c-i) 


(3-1) 


(L-4) 


(T-ll) 


(R-13) 


(H-l) 



Control Systems Using Dynamic Models 



Digital Computer Control 



(A-12) (D-4) (R-20) (L-18) (D- 6 ) 



(S-33) (C-ll) (W-18) 



Hybrid Computer Control 
(D-13) (A-8) 

Analog Computer Control 

(J- 6 ) (L-14) (L- 6 ) (W- 27 ) (H-15) (K-2) (L-3) (R-3) 



Optimal Control 

(L-9) (B-32) (D-13) (A-8) (S-33) (J-5) 

(B-ll) (A-ll) (K-5) (W-29) 

Distributed or Modal Control 

(D-9) (J-6) (M-5) (S-31) (C-4) (S-2) (F-16) (G-13) 
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B2.4 MATHEMATICS AND COMPUTATION 



Ordinary Dlrrerential Equation Theory 



(H-ll) 


(S-21) 


(B-14) 


(1-2) 


(H-12) 


(S-10) 


(H-6) 


Partial Differential Equate 


.on Theory 






(G-7) 


(B-21) 


(0-1) 


(F-16) 


(W-7) 






(F-ll) 


(A-15) 


(F-12) 


(G-6) 


(H-28) 






Integral Equation ' 


Theory 










(T-3) 


(L-8) 


(D-5) 


(s-11) 


(V-4) 


(T-10) 


(W-6) 


(P-4) 


(H-13) 


(m-6) 


(S-14) 


(S-2) 


(G-7) 


(W-7) 


Mathematical Transformations 








(T-4) 


(Z-l) 


(S-31) 


(M-16) 


(R-14) 


(C-8) 




Matrix Mathematics 












(A-4) 


(A-12) 


(H-13) 


(J-7) 


(m- 26 ) 






Numerical 


Solution 


Techniques 








(B-30) 


(F-ll) 


(M-22) 


(P-5) 


(H-25) 


(D-9) 


(w-26) 


(A-15) 


(K-l) 


(D-7) 


(A-10) 


(J-7) 


(D-10) 


(m- 26 ) 


(R-24) 


(L-10) 


(S-13) 


(R-4) 


(L-ll) 


(M-2) 


(J-5) 


(D-8) 


(L-l) 


(C-8) 


(P-8) 


(L-22) 


(N-l) 


(D-15) 


(H-28) 


(L-21) 


(H-19) 


(R-29) 


(L-7) 


(R-17) 


(S-12) 


(F-2) 


(B-32) 


(M-4) 


(P-12) 


(A-17) 


(T-10) 




Boundary Value Problems 










(L-ll) 


(B-8) 


(s-io) 


(H-12) 


(L-7) 


(M-22) 




(H-ll) 


(B-14) 


(V-l) 


(H-25) 


(B-9) 


(P-5) 




(K-l) 


(B-32) 


(F-2) 


(L-22) 


(C-6) 
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Eigen - Values, Vectors, and Functions 



(B-30) 


(H-ll) 


(B-14) 


(M-2) 


(S-31) 


(C-6) 


(0-7) 


(A-12) 


(H-13) 


(0-13) 


(S-13) 


(v-D 


(B-5) 




(G-13) 


(1-2) 


(D-3) 


(S-10) 


(W-26) 


(H-12) 




Special Functions 












(A-18) 


(L-23) 


(m-16) 


(R-28) 


(H-2?) 


(S-25) 


(S-26) 


Computatior 


i and Computer ‘ 


Programming 






(M-4) 


(1-3) 


(0-13) 


(M-18) 


(B-15) 






(0-4) 


(P-8) 


(H-29) 


(s-32) 


(M-25) 
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CHAPTER B3 



BIOGRAPHICAL NOTE 

Lt. Michael N. Hayes was born in Abilene, Texas, on September 2, 
1941, Pour years later the family moved to Ohio, where he attended 
primary schools in Wakeman and in Massillon. He was graduated from 
Washington High School, Massillon, Ohio, in June 1959. 

Mr. Hayes enlisted as a Seaman Recruit in the United States Navy 
in September, 1959* He was graduated from the U.S. Naval Guided 
Missiles School, Dam Neck, Virginia, in July, I960, and was assigned 
to the U.S.S. Kitty Hawk (CVA- 63 ) until June, 1961. He was graduated 
from the U.S. Naval Preparatory School, Bainbridge, Maryland, in 
August, 1961. 

In September, 1961, Mr. Hayes began 4 years of undergraduate 
education at North Carolina State University, Raleigh, N.C., under 
the auspices of the Naval Enlisted Scientific Education Program 
(NESEP), He was graduated in June, 1965, with the degrees of 
Bachelor of Science in both Electrical Engineering and Applied 
Mathematics and with the degree of Master of Electrical Engineering. 

After attending the U.S. Naval Officer Candidate School at 
Newport, R.I., Mr. Hayes was commissioned an Ensign in the U.S, Navy 
in October, 1965, and assigned for duty at the Boston Naval Shipyard, 
Boston, Massachusetts, While stationed at the Shipyard, he worked 
as a Ship Superintendent, supervising naval ship construction and 
repair, and attended evening classes at Boston University in business 
administration and at Northeastern University in pure mathematics. 
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In June, 1967, Mr, Hayes entered M.I.T. as a U.S. Naval Post- 
graduate student in Curriculum 13A with specialized studies in 
electrical engineering. He is a member of Tau Beta Pi, Eta Kappa Nu, 
Phi Eta Sigma, and is an Eagle Scout, 

Lt, Hayes is married to the former Ann Marie Graney of Stoughton, 
Massachusetts. They have a one-year-old son, Stephen, and are 
expecting their second child in December, 1969* 
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